• На главную

            08.02.22 Лекція 1. Метод розв'язання нелінійних рівнянь

            Ссыль на запись: https://youtu.be/33170Mi2gEk

            чисельні розвязання нелінійних рівнянь вигляду

            (1)f(x)=0

            полягає в знаходженні такого x=x* для якого рівняння (1) задовольняється з заданою точністю Е тобто

            (2)|f(x)|E

            Складається з наступних етапів:

            1. відокрремлення кореня
            2. уточнення кореня за допомогою обчислюваного алгоритма

            на першому етапі застосовуєжться графічний або аналітичний підходи

            Аналітичний підхід

            теорема

            неперервна строго монотонна функція f(x) має і при тому єдиний нуль на відрізку [a,b] тоді і тільки тоді коли на його кінцях вона приймає значення різних знаків.

            Графічний спосіб

            побудування графіка функції f(x)

            1. Метод половинного ділення

            попередньо граф методом визначений інтервал локалізації кореня x:x[a0,b0]

            Алгоритм методу:

            1. початковий інтервал [a0,b0] ділиться навпіл: c=a0+b02

            2. Визначаються значення функції f(x) в точках x1=a0, x2=c, x3=b0, тобто f(a0), f(c) f(b0).

            3. Далі обирається та половина початкового інтервалу [a0,b0], на кінцях якої функція f(x) має протилежні знаки. Якщо f(a0)f(c)<0 то наступний відрізок [a1,b1]=[a0,c], інакше [a1,b1]=[c,b0].

            4. Перевіряється умова закінчення процесу пошуку:

              (3)|b1a1|2E

              E - задана точність. Якщо умова (у) виконана, то x=b1+112 , якщо умова (3) не виконана, процес пошуку треба повторити.

            Таким чином на якомусь етапі визначено корінь x* з заданою точністю на основі побудованих вкладених один в одного відрізків[a0,b0],[a1,b1],[a2,b2],..,[ak,bk],bk>ak Для цих відрізків довжини:

            |b1a1|=12|b0a0||b2a2|=12|b1a1|...|b3a3|=12|b2a2|

            звідки і назва методу.

            image-20220208193856276

            рис.1 Метод половинного ділення

            Приклад

            Розв'язати рівняння e2x+3x4=0 з точністю E=102 попередньо графічним способом знайдений інтервал локалізації кореня [a0,b0]=[0,4;0,6]

            k ak bk f(ak) f(bk) ak+bk2
            0 0,4000 0,6000 -0,5745 1,2101 0,5000
            1 0,4000 0,5000 -0,5745 0,2183 0,4500
            2 0,4500 0,5000 -0,1904 0,2183 0,4750
            3 0,4500 0,4750 -0,1904 0,0107 0,4625
            4 0,4625 0,4750 -0,0906 0,0107 0,4688
            5 0,4688 0,4750 -0,0402 0,0107 0,4719
            6 0,4719 0,4750 -0,0148 0,0107 0,4734
            7 0,4734 0,4750 -0,0020 0,0107 [0,4742]

            Табл. 1

            2. Метод Ньютона (метод дотичних)

            Графічна інттерпретація метода image-20220208195255967

            Проведемо дотичну до лінії y=f(x) в точці x=x0. Рівняння дотичної:yy0=f(x0)(xx0) Знайдемо точку перетину дотичної з віссю Oxy0=f(x0)(x1x0) Тобто x1=x0y0f(x0)=x0f(x0)f(x0) Процес повторюється:

            (4)$xk+1=xkf(xk)f(xk);x=0,1,2.

            Для початку обчислень треба задати початкове приближення x0. Умова закінчення процесу пошуку:|x1x0|E Якщо умова (5) виконується, то x=xk+1

            Теорема

            Нехай на відрізку [a0,b0] функція f(x) має першу і другу похідні постійного знаку і нехай f(a0)*f(b0)<0. Тоді, якщо для x0[a0,b0] виконується

            (6)f(x0)f(x0)>0

            то почати з точки x0 послідовність {Xk},(k=0,1,2...), яка задається формулою (4), монотонно збігається до кореня x* , рівняння f(x)=0.

            Приклад

            e2x+3x4=0 з точністю E=102 Маємо:

            f(x)=e2x+3x4f(x)=2e2x+3f(x)=4e2xf(x)таf(x)

            додатні на всій області визначення f(x).

            Виберемо x0=0,6 Маємо:f(0.6)f(0.6)>0 Результат обчислень:

            kxkf(xk)f(xk)f(xk)f(xk)
            0 0,6000 1,1201 9,6402 -0,1162
            1 0,4838 0,0831 8,2633 -0,0101
            2 0,4738 0,0005 8,1585 -0,0001
            3 [0,4737]

            3. Метод січних

            Застосування методу Ньютона потребує на кожній ітерації обчислення значення функції та її похідної. Замінемо похідну функції наближенням різницевим співвідношенням: f(xk)=f(xk)f(xk1)xkxk1 Тоді формула (4) набуває вигляду:

            (7)xk+1=xkf(xk)(xkxk1)f(xk)f(xk1);k=0,1,2..

            Після перетворень:

            (8)xk+1=xk1f(xk)xkf(xk1)f(xk)f(xk1)

            image-20220208201554777

            Легко отримати, що x2=x0y1x1y1y1y0 це точка перетину січної AB з лінією абсцис. Для початку розрахунків за формулою (8) необхідне завдання двох початкових точок x0 x1. Вибір точки х0 можна здійснювати за тим же принципом, що і в методі дотичних. Друга початкова точка х1 обирається в безпосередній близості від х0, бажано між точкою х0 і шуканим коренем.

            Умова закінчення процесу:|xk+1xk|E

            Приклад

            e2x+3x4=0

            x0=0.6;x1=0.59

            k xK f(xk)
            0 0,6000 1,1201
            1 0,5900 1,0244
            2 0,4830 0,0765
            3 0,4744 0,0056
            4 [0,4737]

             

            4. Метод простої ітерації

            Початкове рівняння f(x)=0 приводиться до вигляду

            (9)x=φ(x)

            Розв'язок шукається на основі побудування послідовності

            (10)xk+1=φ(xk);k=0,1,2..

            починаючи з деякого х0.

            Умови збіжності методу і оцінка похибки результату визначаються теоремою:

            Теорема

            Нехай u(x) - визначена і диференційована на [a,b]. тоді, якщо виконується умови 1)φ(x)[a,b],x[a,b] 2)q:|φ(x)q<1|,x(a,b) то рівняння (9)має і при тому єдиний на [a,b] корінь x* і до цього кореня збігається послідовність {Xk},(k=0,1,2...), побудована за формулою (10):limkxk+1=x, яка починається з будь-якого x0[a,b]

            Умови закінчення:|xk+1xk|E

            Зауваження

            Якщо безпосереднє перетворення початкового рівняння f(x)=0 до вигляду х=u(x) не дозволяє отримати умови збіжності методу, можна застосувати таке еквівалентне рівняння:

            (11)x=xλf(x)

            Де λ - параметр, який підбирається таким чином, що в потрібній області виконується умова|φ(x)|=|1λf(x)|q<1

            Приклад

            e2x+3x4=0;x[0.5;0.6]

            Перетворимо початкове рівняння до вигляду

            x=ln(43x)2

            Оскільки для (12) маємо:

            φ(x)=32(43x)maxx[0.5;0.6]32(430.6)=0.6818<1 то процес xk+1=ln(43xk)2,k=0,1,2...збігається. В якості початкової точки візьмемо х0=0,5.

            image-20220208203152947

            15.02.22 Лекція 2. Наближення функції

            Запись: https://youtu.be/sD3HdA475sk

            Наближення функції

            Постановка задачі

            Нехай на деякій множині задана система функцій

            (1)φ0(x),φ1(x),..,φm(x)

            які будемо вважати достатньо гладкими (тобто неперервно диференційованими). Система функцій (1) називається основною системою

            Виведемо поліном

            (2)Qm(x)=c0φ0(x)+c1φ1(x)+..+cmφm(x)

            де C0,C1,C2,...,Cm - постійні коефіцієнти.

            Qm(x) називається узагальненим поліномом порядку m. Наприклад, якщо

            φ0(x)=1;φ1(x)=x;φ2(x)=x2,...,φm(x)=xm, то Qm(x)=c0+c1x+c2x2+...+cmxm маємо звичайний алгебраїчний поліном степені m.

            Якщо

            φ0(x)=1,φ1(x)=cosx,φ2(x)=sinx,φ3(x)=cos2x,φ4(x)=sin2x,Qn(x)=a0+a1cosx+b1sinx+,amcosmx+bmsinmx

            Це тригонометричний поліном порядку m.

            Сформулюємо задачу наближення функції.

            задача

            Задану функцію f(x) треба замінити узагальненим поліномом Qm(x) заданого порядку m так, щоб відхилення функції f(x) від Qm(x) на вказаній множині х було найменшим.

            Q(x) називається апроксимуючим поліномом.

            Множина х може складатися з окремих точок х0,х1,...,хn - тоді маємо точкове неближення.

            Якщо х - це відрізок a<=x<=b, то наближення називається інтегральним

            Відхилення розуміється по-різному: якщо значення f(Xk) i Qm(Xk), k=0,1...,n співпадають - маємо задачу інтерполяції:

            Існують також задачі середньоквадратичного наближення, рівномірного наближення і т.д.

            Інтерполювання функцій

            Нехай у=f(x) задана в точках х0,х1,...,хn своїми значеннями у0,у1,...yn (рис1)

            image-20220216202536153

            рис1

            Будемо вважати функцію f(x) і поліном Qm(x)=c0+c1x+c2x2+...+cmxm близкими, якщо вони співпадають на заданій множині точок х0,х1,...,хn. Ці точки називаються вузлами інтерполяції.

            Задача

            Знайти поліном Qm(x) найнижчої степені m такий, щоQm(xk)=f(xk);k=0,1,2,..

            Qm(x) - називається інтерполяціонним поліномом

            Доведено, що існує єдиний поліном степені не вище n, який приймає в точках х0,х1,...,хn задані значення. Покладемо m=n

            Коефіцієнти полінома Qn(x) знайдемо із системи

            (1){c0+c1x0+c2x02++cnx0n=y0c0+c1x1+c2x12++cnx1n=y1............c0+c1xn+c2xn2+...+cnxnn=yn

            Визначник системи (1) - це визначник Вандельмонда.

            Δ=|1x0x02x0n1x1x12x1n1xnxn2xnn|=(xqxp)0(0p<q<n)

            Система (1) має єдиний роз'вязок. Поліном Ln(x) - це інтерполяційний поліном Лагранжа.

            Він має вигляд:

            (2)Ln(x)=i=0n(xx0)(xx1)(xxi1)(xxi+1)(xxn)(xix0)(xix1)(xixi1)(xixi+1)(xixn)yi

            Отже

            (3)Ln(x)=(xx1)(xx2)(xxn)(x0x1)(x0x2)(x0xn)y0+(xx0)(xx2)(xxn)(x1x0)(x1x2)(x1xn)y1+++(xx0)(xx1)(xxn1)(xnx0)(x11x1)(xnxn1)yn

            Формула (2) - це інтерполяційна формула Лагранжа.

            Приклад

            Знайти поліном степені m<=2, який в точках х0=1, х1=5, х2=-2 приймає значення у0=4, у1=1, у2=7

            Маємо:

            L2(x)=(x5)(x+2)434+(x1)(x+2)471++(x1)(x5)(3)(7)7=13(x5)(x+2)+128(x1)(x+2)++13(x1)(x5)

             

            Приклад 2

            image-20220404220151314

            image-20220404220224936

            Недолік полінома Лагранжа

            при додаванні одного нового вузла треба все перерахувати.

            Роздільні різниці

            Нехай f(x) задана значеннями у0,у1,...,уn в вузлах х0,х1,...хn Складемо відношення:

            [y;x0,x1]=y1y0x2x1

            це число називається роздільною різницею першого порядку функції f(x) в вузлі х0

            Для точки х1 будемо мати: [y;x1,x2]=y2y1x2x1

            Для точки х2: [y;x2,x3]=y3y2x3x2і т.д.

            Роздільна різниця 2го порядку в точці х0:

            [yj;x0,x1,x2]=[yjx1,x2][yjx0,x1]x2x0==y2y1x2x1y1y0x1x0x2x0,

             

            Приклад

            Знайти роздільні різниці для функції у=х^3 в точках х0=1,х1=3,х2=5,х3=2,х4=4

            Складемо таблицю для обчислення різниць

            image-20220216204028034

            Очевидно, що різниці 4го порядку дорівнюють 0 (як і відповідні похідні функції у=х^3)

            Інтерполяційна формула Ньютона

            Запишемо за допомогою роздільних різниць поліном Pn(x), який в точках (вузлах) х0,х1,...,хn приймає значення у0,у1,...,уn Запишемо роздільну різницю [y;x0,x1], а також [y;x0,x1,x2], ..., [y;x0,x1,...,xn] Це ми можемо зробити. Очевидно, що

            [y;x0,x0]=Pn(x0)Pn(x)x0x==Pn(x)Pn(x0)xx0=Pn(x)y0xx0

            Звідки маємо:

            (1)Pn(x)=y0+[y;x,x0](xx0)

            Аналогічно:

            [y;x,x0,x1]=[y;x0,x1][y;x,x0]x1x0=[y;x,x0][y;x0,x1]xx1

            Тобто

            (2)[y;x,x0]=[y;x0,x1]+[y;x,x0,x1](xx1)

            Підставимо (2) в (1). Маємо:

            (3)Pn(x)=y0+[yj;x0,x1](xx0)+[y;x,x0,x1](xx0)(xx1)

            Далі:

            [y;x,x0,x1,x2]=[y;x0,x1,x2][y;x,x0,x1]x2x==[y;x,x0,x1][y;x0,x1,x2]xx2

            Тобто:

            (4)[y;x,x0,x1]=[y;x0,x1,x2]+[y;x,x0,x1,x2](xx2)

            Підставимо (4) в (3):

            (5)Pn(x)=y0+[yjx0,x1](xx0)+[yjx0,x1,x2](xx0)(xx1)++[y;x,x0,x1,x2](xx0)(xx1)(xx2)

            і т.д. В результаті отримаємо:

            (6)Pn(x)=y0+[y;x0,x1](xx0)+[y;;x0,x1,x2](xx0)(xx1)+++[y;x0,x1,,xn](xx0)(xx1)(xxn1)

            Формула (6) - інтерполяційна формула Ньютона. (Відстань між вузлами може бути різна!)

            При появі нового вузла в формулі (6) з'явиться ще один доданок.

            Приклад

            Знайти поліном P3(x), який в вузлах х0=-1,х1=3,х2=4,х3=0 приймає значення у0=4, у1=2 у2=-3, у3=1

            Запишемо обчислення в таблицюimage-20220216212637570

            Маємо:

            (7)P3(x)=412(x+1)910(x+1)(x3)13/30(x+1)(x3)(x4)

            Поліном P3(x) має бути таким же, як і поліном L3(x) Лагранжа.

            Зауваження

            Вузли інтерполяції х0,х1,...,хn можуть бути розташовані в довільній послідовності. Зокрема, якщо вузли перенумеровані в зворотному порядку, поліном (6) приймає вигляд:

            (8)Pn(x)=yn+[xn;xn,xn1](xxn)+[yjxn,xn1,xn2](xxn)(xxn1)++[y;xn,xn1,,x0](xxn)(xxn1)...(xx1)

            Виразу (6) і (8) еквівалентні.

            15.02.22 Лекція 2.5. Кінцеві різниці

            Кінцеві різниці

            Нехай функція y=f(x) задана своїми значеннями в вузлах х0,х1,...,хn. де х1=х0+h, x2=x0+2h, ...

            Тобто вузли є рівновіддаленими. Відповідні значення функції y0, y1, ..., yn

            Визначення

            Кінцевою різницею 1 порядку функції y=f(x) в точці х0 називається величина

            (1)Δf(x0)=f(x1)f(x0),Δy0=y1y0

            Аналогічно вводяться різниці 2го порядку Δ2y0=Δ[Δy0]=Δ[y1y0]=y22y1+y0

            Для точки х1 будемо мати: Δy1=y2y1Δ2y1=y32y2+y1

            Приклад

            Знайти кінцеві різниці для функції у=x^3 в вузлах х0=0, х1=1, х2=2, х3=3, х4=4

            Обчислення зведемо в таблицю

            image-20220216215043189

            Кінцеві різниці image-20220216215147730 називаються кінцевими різницями 1го роду.

            Кінцеві різниці 2го роду позначаються image-20220216215230965 (набла) і обчислюються так:image-20220216215300708

            Тобто береться різниця "назад".

            Різниця 2го порядку:

            image-20220216215331771

            Приклад

            Обчислити кінцеві різниці 2го роду для функції y=x^3 в вузлаї х0=0, х1=1, х2=2, х3=3, х4=4.

            Розрахунки представлені в таблиці

            image-20220216215455366

            Формула Ньютона для рівновіддалених вузлів

            Розглянемо рівновіддалені вузли: х0,х1=х0+h, x2=x0+2h, ... , xn=x0+n*h

            Отримаємо для цього випадку співвідношення між роздільними і кінцевими різницями:[y;,x0,x1]=y1y0x1x0=Δy0h[y;x0,x1,x2]=[y;x1,x2][y;x0,x1]x2x0=

            (3)=y2y1hy1y0h2h=y22y1+y02h2=Δ2y02h2

            Аналогічно неважко отримати, що

            [y;x0,x1,x2,x3]=Δ3y03!h3,,[y;x0,x1,,xn]=Δny0n!hn

            Перепишемо інтерполяційну формулу Ньютона в кінцевих різницях 1го роду:

            (4)Pn(x)=y0+Δy01!h(xx0)+Δ2y02!h(xx0)(xx1)++Δny0n!hn(xx0)(xx1)(xxn1)

            Це інтерполяційна формула Ньютона 1го роду

            Приклад

            Знайти P3(x) за формулою (4), якщо х0=0, х1=2, х2=4, х3=6, у0=1, у1=4, у2=2, у3=7 (Очевидно, що вузли є рівновіддаленими і h=2)

            Розрахунки кінцевих різниць зведемо в таблицю

            image-20220216220005372

            Використовуємо перший рядок (тут більше всього інформації про кінцеві різниці)

            image-20220216220057344image-20220216220109177

            Отримаємо співвідношення між кінцевими різницями 2го роду і відповідними розділеними різницями.

            Маємо:

            image-20220216220213370

            Можна показати, що

            image-20220216220235922

            Отримаємо інтерполяційну формулу Ньютона 2го роду (з використанням кінцевих різниць 2го роду)

            Після підстановки виразів для кінцевих різниць в формулу

            image-20220216220358499

            Отримаємо, що

            Pn(x)=yn+yn1!h(xxn)+yn2!h2(xxn)(xxn1)+++nynn!hn(xxn)(xxn1)(xx1)

            Формули Ньютона (4), (5) - це кінцево-різницевий аналого формули Тейлора.

            Приклад

            Знайти P3(x) за формулою (5), якщо х0=0, х1=2, х2=4, х3=6, у0=1, у1=4, у2=2, у3=7

            Розрахунки зведемо в таблицю

            xkykyk2yk3yk
            01---
            243--
            42-2-5-
            675712

            (таблиця заповнюється знизу)

            Використовуємо останній рядок

            P3(x)=7+51!2(x6)+72!2=(x6)(x4)++123!23(x6)(x4)(x2)=7+52(x6)+78(x6)(x4)++14(x6)(x4)(x2)

            22.02.22 Лекція 2.6. Інтерполяція кубічним сплайном

            Сплайн (spline) -це гнучка лінійка, яка використовується для креслення кривих ліній.

            Дано

            yi=y(xi), xi (i=0,1,2,...,n)

            Треба

            Наблизити y(x) на інтервалі [xi1,xi],(i=1,2,...,n) кубічним поліномом

            (1)

            ui(x)=ai+bi(xxi1)+ci(xxi1)2++di(xxi1)3,i=1,n

            Так щоб

            ui(x)=yi,i=1,n

            і виконати умови неперервності перших і других похідних ui(x) на границях часткових інтервалів[xi1,xi],i=1,n

            Знайдемо невідомі ai,bi,ci,di:

            Для лівого кінця інтервалу [xi1,xi] (за номером це і-й інтервал) з (1) маємо:

            (2)

            yi1=ai,i=1,n

            Для правого кінця

            (3)

            yi=ai+bihi+cihi2+dihi3,i=1,n

            де hi - відстань між вузлами xi1 i xi , i=1,n, тобто hi=xi,i=1,n .

            Рівності (2) і (3) записані з умови співпадання значень функції y(xi) i y(xi1) із значеннями сплайну в вузлах xi і xi1.

            Таким чином, маємо 2n рівнянь (2) і (3) для відшукування 4n невідомих ai,bi,ci,di .

            Необхідно отримати ще 2n умов.

            Для цього будемо вимагати рівність перших і других похідних поліномів (1) у внутрішніх вузлах (їх кількість n1). Це забезпечить гладкість сукупного сплайну.

            Отже маємо:

            ui+1(xi)=ui(xi),                                             (4) ui+1(xi)=ui(xi),i=1,2,,n1 (5)

            Залежності (4) і (5) дадуть ще 2(n1) додаткові умови для визначення ai,bi,ci,di.

            Отримаємо ці умови. Продиференціюємо (1):

            ui(x)=bi+2ci(xxi1)+3di(xxi1)2,i=1,n1, (6)

            ui(x)=2ci+6di(xxi1),i=1,n1                                 (7)

            Вимагаємо неперервності і гладкості у всіх точках, включаючи вузли. Прирівнюючи ліву і праву похідні для внутрішнього вузла xi, маємо:

            bi+2cihi+3dihi2=bi+1,i=1,n1                                        (8)

            ci+3dihi=ci+1,i=1,n1                                                         (9)

            Необхідні ще 2 умови. Їх отримаємо, якщо припустимо, що на кінцях всього інтервалу [x0,xn] другі похідні дорівнюють нулю:

            en(xn)=0,e1(x0)=0                                                                     (10)

            Це відповідає умові незакріплених кінців сплайну (нульова кривизна лінії).

            З (10) маємо:

            12u1(x0)=c1=0                                                     (11)

            12un(xn)=cn+3dnhn=0                                                                 (12)

            Таким чином, отримали систему лінійних рівнянь (2),(3),(8),(9),(11),(12) для визначення 4n невідомих коефіцієнтів ai,bi,ci,di.

            Цю систему можна розв'язати методом Гаусса. Але набагато простіше розв'язати її методом прогонки.

            Для цього попередньо приведемо цю систему до потрібного вигляду.

            З рівняння (2) маємо:

            ai=yi1,i=1,n

            З рівнянь (9) і (12) визначимо:

            di=(ci+1ci)/3hi,i=1,n1                                                     (13)

            dn=cn/3hn                                                                                            (14)

            Якщо покласти cn+1=0, то систему (13),(14) можна записати одним рівнянням:

            di=(ci+1ci)/3hi,i=1,n                                                            (15)

            де cn+1=0

            Підставимо (15) в (3), одночасно виключаючи ai=yi1, отримаємо:

            bi=1hi(yiyi1)13(2сi+сi+1)hi,i=1,n                            (16)

            Виключимо тепер з (8) величини bi та bi+1 за допомогою (16), збільшуючи у другому випадку індекс на одиницю, а величину di на основі (15).

            В результаті маємо систему лінійних рівнянь відносно ci:

            (17)

            {c1=0hi1ci1+2(hi1+hi)ci+hici+1=3(y1yi1hiyi1yi2hi1),i=2,ncn+1=0

            Матриця системи (17) - трьохдіагональна. такі системи розв'язуються спеціальним методом - методом прогонки. Після того, як знайдені коефіцієнти ci , решту коефіцієнтів нескладно знайти з формул (15) і (16)

            Метод прогонки розв'язання систем лінійних рівнянь з трьохдіагональною матрицею

            Представимо друге рівняння системи (17) у вигляді:

            (18)

            uici+1vic1+wici1=Fi,i=2,n

            де позначено:

            ui=hi,vi=2(hi1+hi),wi=hi1,Fi=3(yiyi1hiyi1yi2hi1)

            Коефіцієнти ui,vi,wi і праві частини Fi можна вважати відомими, оскільки їх можна обчислити заздалегідь.

            Покладемо в (18) i=2 і виразимо з рівняння

            u2c3v2c2+w2c1=F2

            Невідомий коефіцієнт c2 через c3 і c1=0:

            (19)

            c2=u2v2c3+w2C1F2v2

            Перепишемо (19) у вигляді

            (20)

            c2=A2c3+B2, де A2=u2v2;B2=w2c1F2v2

            Далі для i=3 з відповідного рівняння (18) маємо

            u3c4v3c3+w3c2=F3

            Підставляючи сюди c2=A2c3+B2, отримаємо

            u3c4v3c3+w3(A2c3+B2)=F3u3c4(v3w3A2)c3+w3B2=F3

            виразимо c3 через c4:

            (21)

            c3=u3v3w3A2c4+w3B2F3v3w3A2

            Представимо (21) у вигляді

            22, 23

            c3=A3c4+B3,   деA3=u3v3w3A2;  B3=w3B2F3v3W3A2

            Для i=4,5,...,n аналогічно отримуємо:

            {C4=A4C5+B4,C5=A5C6+B5,...............Cn=AnCn+1+Bn, дe A4=u4v4w4A3;B4=w4B3F4v4w4A3;A5=u5v5w5A4;B5=w5B4F5v5W5A4,....An=unVnWnAn1;Bn=WnBn1FnVnWnAn1.

            Таким чином, для визначення ci маємо дві системи рівнянь: Перша - відносно невідомих Ai, Bi, (i=2,n):

            24

            {A2=u2v2;B2=F2v2Ak=ukvkWkAk1;Bk=WkBk1FkvkWkAk1,   k=3,n

            Друга - відносно невідомих ci  (i=2,n):

            25

            {C1=0,C2=A2C3+B2,C3=A3C4+B3,........Cn=AnCn+1+Bn,Cn+1=0

            Системи (24)(25) - це системи рекурентних лінійних рівнянь, які розв'язуються послідовно:

            В результаті такої ''прогонки'' знаходяться конфіцієнти ci сплайна (1)

            Алгоритм інтерполювання кубічним сплайном

            1. Обчислити коефіцієнти і праві частини системи рівнянь (18): uici+1vic1+wici1=Fi,i=2,n ui=hi,vi=2(hi1+hi),wi=hi1,Fi=3(yiyi1hiyi1yi2hi1)
            2. Розв'язати систему (18) методом прогонки. Для цього спочатку визначити Ai,Bi  (i=2,n) з системи (24), послідовно позв'язуючи перше, друге,..., останнє рівняння. Далі знайти коефіцієнти ci  (i=2,n) з системи (25), розглядаючи рівнння в зворотній послідовності.
            3. Визначити коефіцієнти ai=yi1,i=1,nbi=1hi(yiyi1)13(2ci+ci+1)hi,i=1,ndi=(ci+1ci)/3hi;i=1,n,

            Приклад

            Есть такая таблица данных

            xi0,01,54,05,27,09,3
            yi1,002,524,163,205,106,20

            Надо получить такую формулу для каждого отрезка, т.е. на x=[0;1,5], на x=[1,5;4] и т.д.

            φi(x)=ai+bi(xxi1)+ci(xxi1)2+di(xxi1)3,i=1,2,,5

            Не сложно догадаться, что i это колво отрезков между заданными иксами, т.е. опять таки i1  это у нас для  x=[0; 1,5], i2 от 1.5 до 4 и т.д.

            Задача сводится к тому, что надо найти эти коефициенты a,b,c,d на разных і:

            по алгоритму определяем коефициенты и правые части системы

            u2=h2=4,01,5=2,5;u3=h3=5,24,0=1,2;u4=h4=7,05,2=1,8u5=h5=9,37,0=2,3;w2=h1=1,50,0;w3=h2=2,5;w4=h3=1,2;w5=h4=1,8;v2=2(1,5+2,5)=8;v3=2(2,5+1,2)=7,4;v4=2(1,2+1,8)=6;v5=2(1,8+2,3)=8,2;F2=3(y2y1h2y1y0h1)=3(4,162,522,52,521,01,5)=1,072;F3=3(y3y2h3y2y1h2)=3(3,204,161,24,162,522,5)=4,368;F4=3(y4y3h3y3y2h2)=3(5,103,201,83,204,161,2)=5,567;F5=3(y5y4h4y4y3h3)=3(6,205,102,35,103,201,8)=1,732;

            Система будет иметь такой вид:

            {1,5c1+8c2+2,5c3=1,0722,5c2+7,4c3+1,2c4=4,3681,2c3+6c4+1,8c5=5,5671,8c4+8,2c5+2,3c6=1,732

            Полученная матрица трёхдиагональная, чтобы решить используем метод прогонки

            A2=u2v2=258=0,1813;B2=w2c1F22v2=0+1,0728=0,1340;A3=u3v3w3A2=1,27,42,5(0,3125)=0,1813;B3=w3B2F3v3w3A2=2,5(0,1340)+4,3687,42,5(0,3125)=0,6093;A4=u4v4w4A3=1,861,2((0,1813)=0,3113;B4=w4B3F4v4w4A3=1,2(0,6093)5,5677,42,5(0,3125)=1,0891;A5=u5v5w5A4=2,38,21,8(0,3113)=0,3011;B5=w5B4F5v5w5A4=1,81,0891+1,73198,21,8(0,3113)=0,4833.

            Теперь можно получить коефициенты сплайна ci:

            c5=A5c6+B5=(0,3011)00,4833=0,4833;c4=A4c5+B4=0,3113(0,4833)+1,0891=1,2396;c3=A3c4+B3=0,18131,23960,6093=0,8340;c2=A2c3+B2=0,3125(0,8340)0,1340=0,1267.

            Знайдемо коефіцієнти ai :

            a1=y0=1,0;a2=y1=2,52;a3=y2=4,16;a4=y3=3,20;a5=y4=5,10.

            Отлично, теперь ищем коефы для b

            b1=1h1(y1y0)13(2c1+c2)h1=11,5(2,521,0)13(20+0,1266)1,5=0,950;b2=1h2(y2y1)13(2c2+c3)h2=12,5(4,162,52)13(20,12660,8340)2,5=1,140;b3=1h3(y3y2)13(2c3+c4)h3=11,2(3,204,16)13(2(0,8340)+1,2396)1,2==0,6286;b4=1h4(y4y3)13(2c4+c5)h4=11,8(5,103,20)13(21,23960,4833)1,8=0,950;b5=1h5(y5y4)13(2c5+c6)h5=12,3(6,205,10)13(2(0,4833)+0)2,3=1,219;с формулок получаем дальшеd1=(c2c1)/3h1=0,1266031,5=0,028;d2=(c3c2)/3h2=0,83400,126732,5=0,1281;d3=(c4c3)/3h3=1,2396+0,834031,2=0,576;d4=(c5c4)/3h4=0,48331,239631,8=0,3191;d5=(c6c5)/3h5=0(0,4833)323=0,070.

            Интерполяция окончена, олды получены

            φ1(x)=1,0+0,950x+0,028x3 для x[0;1,5]φ2(x)=2,52+1,14(x1,5)+0,127(x1,5)20,128(x1,5)3 для x[1,5;4]φ3(x)=4,160,629(x4)0,834(x4)2+0,576(x4)3 для x[4;5,2]φ4(x)=3,200,142(x5,2)+1,240(x5,2)20,319(x5,2)3 для x[5,2;7]φ5(x)=5,10+1,219(x7)0,483(x7)2+0,070(x7)3 для x[7;9,3]

            Тот же пример, но как я понял

             

            0123
            xi1247
            yi2314

            Коефы ищутся исходя из правил

            Вернёмся к формуле

            φi(x)=ai+bi(xxi1)+ci(xxi1)2+di(xxi1)3,i=1,2,,5

            это формулка для каждого отдельного участка сплайна.

            xi1 это начальный х в отрезке сплайна

            теперь попробуем составить это уравнение для первого участка, он же x=[0;1] S0,S1,S2 это названия функций для соответствующих отрезков

            S0=a0+b0(xx0)+c0(xx0)2+d0(xx0)3S1=a1+b1(xx1)+c1(xx1)2+d1(xx1)3S2=a2+b2(xx2)+c2(xx2)2+d2(xx2)3

            начинается этот участок на х=1, знач заменим в формулке x0 на 1:

            S0=a0+b0(x1)+c0(x1)2+d0(x1)3S1=a1+b1(x2)+c1(x2)2+d1(x2)3S2=a2+b2(x4)+c2(x4)2+d2(x4)3

            1. Сплайн должен проходить через узловые точки

            Это значит, что при х=1 у=2, а при х=2 у=3. Из этого ищем коеф a0

            Подставим значения первой и второй точек в формулку, получим:

            S0:2=a0+b0(11)+c0(11)2+d0(11)33=a0+b0(21)+c0(21)2+d0(21)3S1:3=a1+b1(22)+c1(22)2+d1(22)31=a1+b1(42)+c1(42)2+d1(42)3S2:1=a2+b2(44)+c2(44)2+d2(44)34=a2+b2(74)+c2(74)2+d2(74)3

            Это можно упростить

            S0:2=a03=a0+b0+c0+d0S1:3=a11=a1+2b1+4c1+8d1S2:1=a24=a2+3b2+9c2+27d2

            2. В стыках сплайнов всё должно быть гладенько

            На пересечении отрезков должен быть плавный переход, для этого можно использовать производные, так как они показывают динамику роста функции

            Сплайн должен заходить с одинаковой кривизной и углом. Первая производная отвечает за угол. Кривизна же идёт по второй производной, а вот собственно и они:

            S0=(a0+b0(x1)+c0(x1)2+d0(x1)3)=b0+2c0(x1)+3d0(x1)2S1=(a1+b1(x2)+c1(x2)2+d1(x2)3)=b1+2c1(x2)+3d1(x2)2S2=(a2+b2(x4)+c2(x4)2+d2(x4)3)=b2+2c2(x4)+3d2(x4)2(S0=(b0+2c0(x1)+3d0(x1)2)=2c0+6d0(x1)S1=(b1+2c1(x2)+3d1(x2)2)=2c1+6d1(x2)S2=(b2+2c2(x4)+3d2(x4)2)=2c2+6d2(x4)

            В первой точке стыка, она же (2,3) производные сплайнов должны быть одинаковыми, т.е.:

            b0+2c0(21)+3d0(21)2=b1+2c1(22)+3d1(22)22c0+6d0(21)=2c1+6d1(22)

            оно же

            b0+2c0+3d0b1=02c0+6d02c1=0

            аналогично будет в точке (4,1):

            b1+2c1(42)+3d1(42)2=b2+2c2(44)+3d2(44)22c1+6d1(42)=2c2+6d2(44)

            оно же

            b1+4c1+12d1b2=02c1+12d12c2=0

            остнаётся ещё 2 уравнения

            Надо задать поведение сплайнов в начальной и конечной точках

            Зададим нулевую кривизну для начала и конца, т.е.

            S0=0 2c0+6d0(11)=0 2c0=0

            S2=0 2c2+6d2(74)=0 2c2+18d2=0

            Из всего сверху получаем систему из 12 линейных уравнений

            a0=2a0+b0+c0+d0=3a1=3a1+2b1+4c1+8d1=1a2=1a2+3b2+9c2+27d2=4b0+2c0+3d0b1=02c0+6d02c1=0b1+4c1+12d1b2=02c1+12d12c2=02c0=02c2+18d2=0

            учитывая, что мы знаем некоторые коэфы - заменяем известные:

            a0=2a1=3a2=1c0=0b0+d0=12b1+4c1+8d1=23b2+9c2+27d2=3b0+3d0b1=06d02c1=0b1+4c1+12d1b2=02c1+12d12c2=02c2+18d2=0

            Таким образом получаем матрицу СЛАУ значений коефов А и матрицу столбец значений В

            30.03.22 Лекція 3. Середньоквадратичне наближення функцій

            Раніше було розглянуто наближення функцій методами інтерполювання. Вони забезпечують хорошу точність на невеликих відрізках. Середньоквадратична апроксимізація дозволяє будувати наближені формули у випадку великих відрізків з великою кількістю вузлів

            Задача середньоквадратичного наближення

            Для заданої функції f(x) треба знайти поліном

            (1)Qm(X)=c0μ0(x)+c1μ1(x)+...+cmμm(x)

            де μ0(X),μ1(x),..,μm(x) - задана система функцій, таких щоб середньоквадратичне відхилення від заданої функції f(х) було мінімальним. Середньоквадратичним відхиленням називається усереднення на заданій множині точок квадрата різниці заданої функції f(x) і полінома Qm(X). У випадку, коли функція f(x) задана своїми значеннями у вузлах x0,x1,..,xn за міру близості можна взяти величину E, яка виражається формулою:

            (2)E=i=0n(f(xi)Qm(xi))2

            Якщо множина Х, на якій розглядається задача середньоквадратичного наближення, задана у вигляді інтервалу X=(x|axb),

            (3)E=ab(f(x)Qm(x))2dx

            Принципова відмінність задачі середньоквадратичного наближення від задачі інтерполяції полягає в тому, що кількість невідомих коефіцієнтів c0,c1,..,cm суттєво менше, ніж кількість вузлів апроксимації, тобто m<<n.

            image-20220330193638128

            Якщо розглядається міра відхилення 2, маємо задачу точкової середньоквадратичної апроксимації, якщо має місце міра відхилення 3 - це неперервна апроксимація.

            Точкове середньоквадратичне наближення функцій

            В даному випадку задача розглядається для функції f(x), яка задана своїми значеннями y0,y1,...,yn в вузлах x0,x1,..,xn, а міра відхилення - формулою 2.

            Оскільки E є функцією невідомих коефіцієнтів c0,c1,..,cm, умова мінімізації оцінки 2 має вигляд

            (4)Eci=0,i=0,1,,m

            Розглянемо випадок алгебраїчного поліному

            Qm(x)=c0+c1x+c2x2++cmxm. тоді

            E=i=0n(c0+c1xi+c2xi2++cmximyi)2.

            Візьмемо похідні 4, отримаємо в результаті систему (m+1) алгебраїчних рівнянь відносно невідомих c0,c1,..,cm

            (5){i=0n[(c0+c1xi+c2xi2++cmx1myi)1]=0i=0n[(c0+c1xi+c2xi2++cmx1myi)xi]=0.......i=0n[(c0+c1xi+c2x12++cmximyi)xim]=0

            система 5 може бути переписана у вигляді

            (6){c0i=0n1+c1i=0nxi+c2i=0nxi2++cmi=0nxim=i=0nyic0i=0nxi+c1i=0nxi2+c2i=0nxi3++cmi=0nxim+1=i=0nyixi...........c0i=0nxim+c1i=0nxim+1+c2i=0nxim+2++cmi=0nxi2m=i=0nyixim

            Матриця системи 6 симетрична відносно головної діагоналі

            (i=0n1i=0nxii=0nxi2i=0nximi=0nxii=0nxi2i=0nxi3i=0nxim+1i=0nximi=0nxim+1i=0nxim+2i=0nxi2m)

            Симетричність матриці є відмінною властивістю застосування методу найменших квадратів, який полягається в основу середньоквадратичного наближення. Можна довести, що якщо серед вузлів x0,x1,..,xn нема співпадаючих, то визначник системи 6 відмінний від нуля, тобто існує єдиний розв'язок системи.

            Якщо m-n, то апроксимуючий поліном Qm(x) співпадає з поліномом Лагранжа для системи вузлів x0,x1,..,xm, причому E=0.

            Таким чином, наближення у вигляді апроксимації є більш узагальненим процесом в порівнянні з інтерполяцією.

            Система 6 називається нормальною системою

            Приклад точкової середньоквадратичної апроксимації

            Нехай функція y=f(x) задана своїми значеннями в вузлах:

            xi0,781,562,343,123,81
            yi2,501,201,122,254,28

            Наблизимо функцію поліномом Q2(x)

            Q2(x)=c0+c1x+c2x2

            Нормальна система рівнянь 6 для знаходження невідомих c0,c1,c2 має вигляд:

            {c0i=041+c1i=04xi+c2i=04xi2=i=04yic0i=04xi+c1i=04xi2+c2i=04xi3=i=04yixic0i=04xi2+c1i=04xi3+c2i=04xi4=i=04yixi2

            Обчислення коефіцієнтів нормальної системи запишемо в табличному вигляді

            image-20220330202039486

            Таким чином маємо систему:

            {5c0+11,61c1+32,768c2=11,35011,61c0+32,768c1+102,761c2=29,77032,768c0+102,761c1+341,750c2=94,604

            Розв'язок цієї системи: c0=5,045,c1=4,043,c2=1,009

            Апроксимуючий поліном має вигляд

            Qz(x)=5,0454,043x1+1,009x2

            Знайдемо значення величини відхилення image-20220330202315842

            ε=2,29104=0,015

            Для алгебраїчного полінома основна система функцій

            φ0(x)=1,φ1(x)=x,φ2(x)=x2,,lnmmi(x)=x3,

            Якщо ввести позначення

            (φ0,φ0)=i=0nε0(εi)φ0(xi)=i=0n1(φk,φj)=i=0nφk(xi)φj(xi)=i=0nxikxi(φk,y)=i=0nφk(xi)y(xi)=i=0nxikyi,

            То систему нормальних рівнянь 6 можна переписати в компактному вигляді

            (7){c0(φ0,φ0)+c1(φ0,φ1)+c2(φ0,φ2)++cm(φ0,φm)=(φ0,y)c0(φ1,φ0)+c1(φ1,φ1)+c2(φ1,φ2)++cm(φ1,φm)=(φ1,y)..........c0(φm,φ0)+c1(φm,φ1)+c2(φm,φ2)++cm(φm,φm)=(φm,y)

            Приклад 2

            Функцію f(x)=1x2π2, задану своїми значеннями в вузлах x0=0,x1=π/3,x2=π/2,x3=2π/3,x4=pi наблизити тригонометричним поліномом

            Q1(x)=C0+C1cosx

            В цьому випадку φ0(x)=1,φ1(x)=cosx.

            Нормальна система має вигляд:

            (8){c0(φ0,φ0)+c1(φ0,φ1)=(φ0,y)c0(φ1,φ0)+c1(φ1,φ1)=(φ1,y)

            Обчислення коефіцієнтів системи 8 запишемо в таблицю:

            image-20220330203502570

            Нормальна система розпадається на два незалежних рівняння

            {5c0=115/36c0=23/362,5c1=7/6c1=7/15

            Тобто Q1(x)=23/36+7/15

            Лекція 5. Основні теореми квадратичного наближення

            Нехай функція f(X) задана на множині x=(x|axb) і вибрана система лінійно незалежних функцій μ0(x),μ1(x),..,μm(x) Необхідно так підібрати коефіцієнти c0,c1,..,cm, щоб поліном

            (1)h(x)=c0μ0(x)+c1μ1(x)+...+cmμm(x)

            був якомога близьким до f(x) на множині вузлів axib в сенсі мінімуму оцінки E

            Розглянемо систему нормальних рівнянь в загальному випадку:

            (1){(φ0,φ0)c0+(φ0,φ1)c1++(φ0,φm)cm=(φ0,f)(φ1,φ0)c0+(φ1,φ1)c1++(φ1,φm)φm=(φ1,f)......(φm,φ0)c0+(φm,φ1)c1++(φm,φm)φm=(φm,f)

            Теорема 1.

            Для того, щоб поліном h(x) був розв'язком задачі найкращого квадратичного наближення, необхідно і достатньо, щоб числа c0,c1,..,cm були розв'язками системи 1.

            Достатність

            Нехай числа c0,c1,..,cm - розв'язки системи 1. Візьмемо деякі інші числа b0,b1,..,bm і складемо функцію

            (2)g(x)=b0φ0(x)+b1φ1(x)++bmφm(x)

            Доведемо, що fg>fh, де позначено f=f(x),f=g(x),h=h(x),

            Маємо:

            (*)fg2=(fg,fg)=(fh+hg,fh+hg)==((fh)+(hg),(fh)+(hg))=fh2++2(fh,hg)+hg2

            Перепишемо перше рівняння системи у вигляді

            (3)(f(c0φ0+c1φ1++cmφm),φ0)=0,де φk=φk(x)

            Приймаючи до уваги 1, отримаємо

            (5){(fh,φ0)=0(fh,φ1)=0........(fh,φm)=0

            умова 5 - це умова ортогональної функції (fh) і φ0,φ1,...,φm.

            Запишемо

            (6)h(x)g(x)=i=0m(eibi)φi(x)

            Сформуємо:

            (fh,hg)=(fh,i=0m(cibi)φi)==(c0b0)(fh,φ0)+(c1b1)(fh,φ1)++(cmbm)(fh)=0,

            тобто

            (7)(fh,hg)=0

            з урахуванням 7 перепишемо рівність * у вигляді

            (8)fg2=fh2+hg2

            Рівняння 8 - це теорема Піфагора в Гільбертовому просторі. Оскільки hg2>0 в умовах hg, то

            fg2>fh2

            що і потрібно було довести.

            Доведемо необхідність.

            Нехай fg>fh. при будь-якому g=g(x). необхідно довести, що числа c0,c1,..,cm задовольняють системі 1

            Запишемо *:

            (8)fg2=fh2+2(fh,hg)+hg2

            Оскільки 8 виконується при будь-якому g=g(x), розглянемо випадок g(x), близького до h(x), тобто hg

            Тоді величина hg2 - є величина другого порядку малості у зрівнянні з величиною 2(fh,hg)

            Тоді в першому наближенні

            (9)fg2=fh2+2(fh,hg)

            Оскільки(fh,hg) може приймати різні знаки, то для виконання 9 має бути

            (10)(fh,hg)=0

            Таким чином з 9 будемо мати

            ||fg||2=||fh||2

            тобто з 10 отримаємо з урахуванням 6 що

            (fh,hg)=(fh,i=0m(cibi)φi(x))==(с0b0)(fh,φ0)+(c1b1)(fh,φ1)++(cmbm)(fh,φm)=0

            Оскільки функція g(x) - довільна, то b0,b1,..,bm - теж довільні, тоді довільні і коефіцієнти (c0b0),(c1b1),..,(cmbm) Тоді з 10 випливає що одночасно

            (11){(fh,φ0)=0(fh,φ1)=0(fh,φm)=0

            А система 11 як ми бачили раніше - це фактично система 1. тобто необхідність доведення.

            Теорема 2.

            Система нормальних рівнянь 1 для будь-якої функції f(x) має розв'язок і при тому єдиний.

            Доведемо від противного. Нехай система 1 має два розв'язки: c0,c1,...,cm і b0,b1,...,bm Тоді маємо:

            (12)h(x)=c0φ0(x)+c1φ1(x)++cnφn(x)

            в скороченому вигляді

            h(x)=i=0mсiφi(x).

            Аналогічно:

            (13)g(x)=b0φ0(x)+b1φi(x)++bmφm(x)==i=0mbiφi(x)

            Оскільки ci,i=1,2,..,m;bi,i=0,1,2,..,m - розв'язки, то одночасно

            (fh,φi)=0 ,i=0,1,..,m (14)

            (fg,φi)=0 ,i=0,1,..,m (15)

            Віднімемо 14 від 15, маємо:

            (fg,φi)(fh,φi)=0,i=0,m(fgf+h,φi)=0,i=0,m.

            Остаточно отримали:

            (16)(hg,φi)=0i=0,m

            Розпишемо 16 для i=0,1,..m:

            (17){(hg,φ0)=0(hg,φ1)=0(hg,φm)=0

            Домножимо перше рівняння 17 на (c0b0), друге - на(c1b1) і т.д. і складемо:

            (c0b0)(hg,φ0)+(c1b1)(hg,φ1)++(cmbm)(hg,φm)=0(hg,(c0b0)φ0+(c1b1)φ1++(cmbm)φm)=0

            тобто (hg,hg)=0, це означає, що ||hg||2=0, звідки випливає, що h=g, тобто ci=bi,i=0,m, що є протиріччям. Теорема доведена

            Похибка квадратичного наближення

            Розглянемо

            (18)fh2=(fh,fh)=f22(f,h)+h2

            Оскільки

            (19)(f,h)=(fh+h,h)=(fh,h)+h2

            Але оскільки (fh,h)=0, бо (fh,φi)=0,i=0,m, то з 18 з урахуванням 19 маємо

            fh2=f22h2+h2=f2h2

            Остаточно для похибки маємо:

            (20)fh2=f2h2

            з 20 випливає що значення похибки заздалегідь невідоме, тобто оцінка 20 є постеріорною (після досліду).

            Неперервне квадратичне наближення алгебраїчними поліномами

            Нехай на множині X=(x|axb) задана система базисних функцій

            φ0(x)=1,φ1(x)=x,φ2(x)=x2,,φm(x)=xm

            Апроксимуючий поліном має вигляд

            hm(x)=c0+c1x+c2x2+1.+cmxm

            Система нормальних рівнянь записується у вигляді

            {(1,1)c0+(x,1)c1++(xm,1)cm=(f,1)(1,x)c0+(x,x)c1++(xm,x)cm=(f,x)............(1,xm)c0+(x,xm)c1++(xm,xm)cm=(f,xm)

            де

            (1,1)=ab11dx=ba(x,1)=ab1xdx=abxdx=x22|ab=12(b2a2)(xm,1)=ab1xmdx=xm+1m+1|ab=bm+1am+1m+1

            Приклад

            задану функцію f(x)=x+1 наблизити поліномом Q(x)=c0+c1x на відрізку xє[0,3]

            запишемо систему нормальних рівнянь:

            (21){03dxc0+03xdxc1=03x+11dx03xdxc0+03x2dxc1=03x+1xdx

            Обчислимо інтеграли в системі 21

            03dx=3,03xdx=x22|03=9203x+1dx=03(x+1)1/2dx=(x+1)3/23/2|03==23(43/21)=14303x2dx=x33|03=273=903x+1xdx=|u=xdu=dxdv=x+1dxv=x+1dx=∣=|==x23(x+1)3/2|032303(x+1)3/2dx=(x+1)3/23/2=32343/202325(x+1)5/2|03==223415(45/21)=16415(321)=1643115=11615

             

            Нормальна система:

            {3c0+92c1=14/392c0+9c1=11615

            Домножимо друге рівняння на 12 і додамо до першого:

            Тоді з другого рівняння отримаємо:

            Таким чином: Q(x)=1615+44135x

            Оцінимо похибку наближення:

            Похибка достатньо мала

            05.04.22 Лекція 6

            Неперервне квадратичне наближення у випадку ортогонального базису

            нехай базисні функції φ0(x),φ1(x),...,φm(x) взаємно ортогональні на [a.b] тобто (φk,φp)=0,kp тоді система нормальних рівнянь має вигляд:

            (1){(φ0,φ0)c0=(f,φ0)(φ1,φ1)c1=(f,φ1).....(φm,φm)cm=(f,φm)

            і розпадається на m окремих рівнянь

            (2)(φk,φk)сk=(f,φk),k=0,m

            з рівнянь 2 маємо:

            (3)ck=(f,φk)(φk,φk),k=0,m

            неперервне квадратичне наближення тригонометричними поліномами

            Розглянемо систему функцій

            φ0(x)=1φ1(x)=cos(x)φ2(x)=sin(x)φ3(x)=cos(2x)φ4(x)=sin(2x)....φ2m1(x)=cos(mx)φ2,(x)=sin(mx)

            ці функції є ортогональними на інтервалі x[π,π] оскільки

            ππcoskxcospxdx=20πcoskxcospxdx==0π(cos(k+p)x+cos(kp)x)dx==sin(k+p)xk+p|0π+sin(kp)xkp|0π=0
            ππsinkxsinpxdx=20πsinkxsinpxdx==0π(cos(kp)xcos(k+p)x)dx=sin(kp)xkp|0πsin(k+p)xk+p|0π=0ππsinkxcospxdx=12ππ(sin(k+p)x+sin(kp)x)dx==12(cos(k+p)xk+p|ππcos(kp)xkp|ππ)==12(cos(k+p)π+cos(k+p)πk+p+cos(kp)π+cos(kp)πkp)==0

            Тригонометричний поліном степені m має вигляд:

            (4)Qm(x)=a0+a1cos(x)+b1sin(x)+a2cos(2x)+b2sin(2x)+...+amcos(mx)+bmsin(mx)
            (1,1)=ππdx=x|ππ=2π(coskx,coskx)=ππcos2kxdx=12ππ(1+cos2kx)dx==12(x+12sin2x)|ππ=π(sinkx,sinkx)=ππsin2xdx=12ππ(1cos2x)dx==12(x12sin2x)|ππ=π

            Система нормальних рівнянь має вигляд:

            {(1,1)a0=(f,1)(coskx,coskx)ak=(f,coskx),k=1,m(sinkx,sinkx)bk=(f,sinkx),k=1,m

            звідки:

            (5)a0=12πππf(x)dx
            (6)ak=1πππf(x)coskxdx
            (7)bk=1πππf(x)sinkxdxk=1,m

            коефіцієнти a0,ak,bk співпадають з коефіцієнтами ряда фур'є для фунції f(x):

            (8)f(x)=a¯02+n=1(ancosnx+bnsinnx), де a¯0=1πππf(x)dx.

            таким чином, ряд фур'є є оптимальним наближення функції на інтервалі x[π,π] в сенсі мінімальності середньоквадратичної похибки.

            Приклад

            функцію f(x)=|x| на інтервалі x[π,π] наблизити тригонометричним поліномом 3го порядку

            Q3(x)=a0+a1cosx+b1sinx+a2cos2x+b2sin2x++a3cos3x+b3sin3x

            знайдемо коефіцієнти:

            a0=12πππ|x|dx=1π0πxdx=1πx22|0π=π2

            отримали: Q3(x)=π24πcosx49πcos3x

             

            оцінимо похибку наближення

            Лекція 7 наближене обчислення похідних і інтегралів

            табличне (чисельне) диференціювання

            Спочатку розглянемо приклад. Є функція f(x)=Esin(xE2) де E>0, яка задовольняє умові |f(x)|E отже є обмеженою

            image-20220405174855092

            Знайдемо похідну f(x)=1Ecos(xE2) отже 1Ef(x)1E. Якщо взяти дуже мале E (близько 0) то для обмеженої функції f(x) похідна буде приймати дуже великі значення.

            Оскільки чисельні методи визначення похідних є наближенними, то застосовувати їх треба з обмеженістю, бо похибка може бути занадто великою. Якщо значення функції задані як результати вимірювань і необхідно дослідити поведінку похідних, то треба застосувати попереднє зглажування з подальшим диференціюванням.

            Формули Грегорі-Ньютона

            Нехай функція f(x) задана своїми значеннями в рівновіддалених вузлах x0,x1=x0+h,x2=x0+2h,.. . Необхідно обчислити наближене значення f(x0) Використаємо інтерполяційну формулу Ньютона: (1)

            Введемо позначення yk=f(xk), xx0h=u, тоді

            (2)

            запишемо з урахуванням (2)

            і т.д. Тоді з 1 маємо:

            (3)

            Групуючи коефіцієнти по степеням та в 3 отримаємо:

             

            Покладемо hu=τ і перепишемо 4 у вигляді:

            (5)

            Порівняємо формулу 5 з розкладанням функції f(x) в ряд Тейлора:

            (6)f(x0+τ)=y0+y01!τ+y02!τ2+...

             

            Оскільки розкладання в ряд по степеням τ можна виконати єдиним способом, то порівнюючи коефіцієнти при однакових степенях τ в правих частинах рівностей 5 і 6 отримаємо:

            (7)y0=(1/h)(Δy0Δ2y02+Δ3y03...)
            (8)y0=(1/h2)(Δ2y0Δ3y0+1112Δ4y0...)

             

            Формули 7 і 8 це формули Грегорі-Ньютона, які використовують кінцеві різниці 1го роду. Вони апроксимують першу і другу похідні функції f(x)

            Формули 7 і 8 записані для вузла x0​. Для довільного вузла xk​ ці формули мають вигляд:

            (9;10)f(xk)=1h(ΔykΔ2yk2+Δ3yk3)f(xk)=1h2(Δ2ykΔ3yk+1112Δ4yk)

             

            Приклад

            Знайти f(2) якщо

            image-20220405180803992

            Знайдемо кінцеві різниці. Результат представимо в таблиці

            image-20220405180836135

            Маємо:

            image-20220405180917665

            Для другої похідної маємо:

            image-20220405180945756

            Лекція 8 13.04.22

            Операторно-символічне виведення формули Грегорі-Ньютона

            Запишемо формулу Тейлора для функції y=f(x)​:

            (11)f(x0+h)=f(x0)+f(x0)1!h+f(x0)2!h2++f(n)(x0)n!hn

             

            Введемо оператор диференціювання:

            (12)Df(x)=f(x)

             

            Тоді f(x)=Df(x)=D2f(x),f(n)(x)=Dnf(x)

            Перепишемо формулу 11 з урахування введеного оператора диференціювання: f(x0+h)=f(x0)+Df(x0)1!h+D2f(x0)2!h2++Dnf(x0)n!hn

            Використовуючи формальні операції отримаємо:

            (13)f(x0+h)=f(x0)(1+D1!h+D22!h2++Dnn!hn)

            Порівняємо вираз в дужках в правій частині 13 з розкладенням функції ex в степеневий ряд:

            Отримаємо:

            (14)f(x0+h)=eDhf(x0).

             

            Віднімемо f(x0)​ з лівої і правої частини 14:

            (15)f(x0+h)f(x0)=eDhf(x0)f(x0)=(eDh1)f(x0).Δf(x0)=(eDh1)f(x0)

             

            Рівність 15 треба розуміти як рівність дій, а не величини. Отже з (15) маємо:

            Δ=eDh1;eDh=1+Δ

            Логарифмуємо:

            (16)Dh=ln(1+Δ),D=1hln(1+Δ).

            Оскільки для ln(1+Δ) маємо розкладання:

            ln(1+Δ)=ΔΔ22+Δ33Δ44+

            то 3 (16)

            =1h(ΔΔ22+Δ33Δ44+) (17) 

            Домножимо (17) на f(x0)б маємо:

            Df(x0)=1h(ΔΔ22+Δ33Δ44+)f(x0)

            тобто

            f(x0)=1h(Δf(x0)Δ2f(x0)2+Δ3f(x0)3Δ4f(x0)4+)

            Отримаємо формулу для f(x0) Возведемо (17) в квадрат:

            D2=1h2(ΔΔ22+Δ33Δ44+)2==1h2(Δ2+Δ44Δ3+23Δ4)==1h2(Δ2Δ3+1112Δ4) Тоді Δ2f(x0)=1h2(Δ2Δ3+1112Δ4)f(x0)

            і остаточно:

            f(x0)=1h2(Δ2f(x0)Δ3f(x0)+1112Δ4f(x0))

            Виведемо формули для обчислення похідних з використанням різниць другого роду. Запишемо формулу Тейлора у вигляді

            f(xnh)=f(xn)f(xn)1!h+f(xn)22!h2++(1)nf(n)(xn)n!hn

            Використаємо оператор диференціювання D:

            f(xnh)=f(xn)f(xn)1!n+2f(xn)2!h2++(1)nDnf(xn)n!hn+f(xnh)=(1Dh1!+D2h22!+(1)nDnhnn!+)f(xn)(18)

            Оскільки

            1ϕh1!+ϕ2h22!=eh,

            то формула (18) набуває вигляду:

            f(xnh)=ehf(xn)

            тоді

            f(xnh)f(xn)=(eDn1)f(xn)

            тобто

            f(xn)=(eh1)f(xn). (19) 

            Прирівнюючи дії в виразі (19) маємо

            eDh1=,eDh=1,

            Логарифмуємо:

            h=ln(1)=223344.

            Отримаємо:

            D=1n(223344)==1n(+22+33+44+)(20)

            Домножимо (20) на f(xn), отримаємо:

             D f(xx4)=1n(+22+33+44+)f(xn)

            тобто

            f(xn)=1h(f(xn)+2f(xn)2+3f(xn)3+)(21).

            Для другої похідної отримаємо:

            f(xn)=1h2(2f(xn)+3f(xn)+11124f(x11)+) (22) 

            Приклад

            Для функції y=x3 отримати значення похідних y(5) та y(5), використовуючи формули (21) та (22).

            Сформуємо таблицю

            xkykyk2yk3yk4ykykyk
            00      
            111     
            2876    
            32719126   
            464371860  
            5125612460  
            y(5)=11(61+242+63=61+12+2=75y(5)=112(24+6)=30

            Оскільки 4yk=0, то формули 21 і 22 дають точні значення похідних y(5) і y(5).

            Наближення обчислення первісних функцій. Перша та друга формули Адамса

            Постановка задачі

            Функція y=f(x) задана в рівновіддалених вузлах x0,x1,x2,..,xn. де xi=x0+ih, своїми значеннями y0,y1,..,yn

            Введемо позначення:

            F(xn+1)=x0xn+1f(x)dx=x0xnf(x)dx++xnxn+1f(x)dx=F(xn)+xnxn+1xn+1f(x)dx.
             Отже F(xn+1)=F(xn)+xnxnxn+1(x)dx

            Квадратичною формулою називається всяка проста формула, що апроксимує інтеграл enxint f(x)dx Замінимо підінтегральну функцію f(x) інтерполяційним многочленом Ньютона P3(x) :

            xnxn+1f(x)dxxnxn+1P3(x)dx=xnxn+1(yn+yn11!h(xxn1)++2yn2!h2(xxn1)(xxn1)+3yn13!h3(xxn1)(xxn1)(xxn2))dx

            Позначимо

            xxnh=u,x=xnH=0x=xn+1u=1

            Тоді:

            xxn=hu.dx=hdu

            Маємо:

            Тоді:

             

            Отримали першу формулу Адамса:

            (2)xnxn+1f(x)dx=h(yn+12yn+512yn2+38yn)

            Це - наближена формула для обчислення інтеграла на інтервалі x[xn,xn+1]. Вона використовує значення функція f(x) за межами цього інтегралу. Це - недолік формули. Отримаємо другу формулу Адамса: для цього використаємо значення функції f(x) у вузлах xn1,xn,xn+1: image-20220413104454444 Введемо заміну image-20220413104512580 тоді: image-20220413104527456 image-20220413104539453

            Друга формула Адамса

            Остаточно:

            (3)xnxn+1f(x)dx=h(yn+112yn+11122yn+11243yn+1)

            Квадратична формула трапецій

            Замінимо підінтегральну функцію f(x) на інтервалі [xk,xk+1] поліномом Ньютона першого порядку (рис1.) Маємо:

            image-20220413104742619

            image-20220413104756930

            image-20220413104808398 Формула (4) - квадратична формула трапецій. На рис 1. зображена трапеція, що апроксимує площину під графіком функції f(x) на інтервалі [xk,xk+1]. Отримаємо складову формули трапецій, що має міск при обчисленні інтеграла image-20220413104940358 Розіб'ємо інтервал [a,b] на n рівних інтервалів і застосуємо до кожного з них квадратичну формулу (4): image-20220413105043395 Отримали велику (складову) формулу трапецій: image-20220413105121628

            Квадратична формула середніх (прямокутників)

            Формулу ілюструє рис 2. image-20220413105213975 Площина під графіком функції y=f(x) замінюється площиною прямокутника зі сторонами h і y(xk+1/2), де xk+1/2=xk+1+xk2. Таким чином: image-20220413105423156 Легко отримати і велику формулу прямокутників image-20220413105513366

            Квадратична формула Сімпсона

            Наблизимо підінтегральну функцію f(x) на інтервалі [xk,xk+1] поліномом 2го порядку. Для цього використаємо вузли xk,xk+1/2=xk+h/2

            image-20220413105701771

            image-20220413105715129

            image-20220413105722070

            Отримали квадратичну формулу Сімпсона. Виведемо велику (складову) формулу Сімпсона. Інтервал [a,b] розіб'ємо на n рівних інтервалів довжиною ban, застосовуючи до кожного з них формулу (8) , маємо

            image-20220413105913236

            image-20220413105920553

             

            21.04.22 Похибки квадратичних формул

            На попередній лекції ми отримали наступні квадратичні формули:

            Формула трапецій:

            image-20220423171341047формула середніх:

            image-20220423171501312

            Формула Сімпсона:

            image-20220423171519327

            Оцінимо похибки квадратичних формул. Розкладемо підінтегральну функцію y=f(x) в ряд Тейлора з центром розкладання x=xk+1/2:

            image-20220423171648142

            Проінтегруємо (1) на інтервалі [xk,xk+1]: image-20220423171736558 image-20220423171828599 Порівнюючи вираз (5) з формулою середніх (2), отримаємо, що похибка квадратичної формули середніх: image-20220423171915702 Щоб отримати похибку формули трапецій, підставимо в формулу (4) x=xk, отримаємо: image-20220423171959094 image-20220423172014558 Якщо підставити в (4) x=xk+1, отримаємо: image-20220423172043762 Враховуючи (7) і (8), маємо: image-20220423172117400 image-20220423172128790

            Віднімаючи (9) від (5), отримаємо похибку квадратичної формули трапецій: image-20220423173000750

            image-20220423173011731

            Таким чином формула середніх вдвічі точніша за формулу трапецій. Причому значення похибки протилежні.

            Розглянемо похибку квадратичної формули Сімпсона. Сформуємо формулу Сімпсона з використанням формул (7) і (8): Маємо:

            image-20220423173244070

            Порівнюючи (11) з розкладанням (5), отримаємо, що image-20220423173432349

            Чисельні та чисельно аналітичні методи розв'язання задачі Коші для звичайних диференціальних рівнянь

            Задача Коші. Загальні відомості

            Розв'язання звичайних диференціальних рівнянь займає важливе місце серед прикладних задач фізики, хімії та техніки.

            Відомо, що звичайне диференціальне рівянння р-го порядку можна за допомогою заміни звести до еквівалентної системи р-рівнянь 1го порядку

            приклад

            image-20220423184656518 image-20220423184707313

            Аналогічно, довільну систему диференціальних рівнянь любого порядку можна замінити деякою еквівалентною системою рівнянь 1го порядку. Тому в подальшому будемо розглядати систему рівнянь 1го порядку

            image-20220423184836004 Представляючи її в векторному вигляді

            image-20220423184915081 Як відомо з теорії диференціальних рівнянь, система (1) має множину розв'язків, яка в загальному випадку залежить від р параметрів c1,c2,..,cp. Введемо вектор image-20220423185039613 Тоді розв'язок системи (1) або векторного рівняння (2) можна записати у формі: image-20220423185113562 Для визначенння c1,c2,..,cp, тобто для виділення єдиного розв'язку, необхідно накласти р додаткових умов на функції yk(x).

            Розрізняють 3 типи задач для звичайних диференціальних рівнянь:

            1. Задачі Коші
            2. крайтова задача
            3. задача на власні значення

            Задача Коші (задача з початковими умовами)

            Маємо додаткові умови для визначення постійних c1,c2,..,cp: image-20220423185354157 (тобто задані значення всіх функцій yk(x) в одній і тій же точці x=x0).

            Ці умови можна розглядати як завдання координат початкової точки (x0,y0,y20,...,yp0) інтегральної кривої в (p11)-мірному просторі (x,y1,y2,..,yp) Розв'язок рівняння (2) при цьому необхідно знайти на деякому відрізку x[x0,xn]

            Треба зауважити, що якщо праві частини рівнянь (1) неперервні і обмежені в деякому околі початкової точки (x0,y10,y20,..,yp0), то задача Коші image-20220423193606195 Має розв'язок, але взагалі кажучи, не єдиний.

            Якщо праві частини f1,f2,..fp не тільки неперервні, але й задовольняють умові Ліпшиця по змінним yk, то розв'язок задачі Коші єдиний.

            Умова Ліпшиця:

            image-20220423193827936

            Методи розв'язання задачі Коші

            Поділяють на точні, наближені та чисельні. (точні методи розглядаються в курсі диференціальних рівнянь)

            Наближені методи

            Розв'язок отримується як границя y(x) деякої послідовності yn(x), причому yn(x) виражаються через елементарні функції або за допомогою квадратур. Обмежуючись кінцевим числом n, отримаємо наближений вираз для y(x). Приклад: метод розкладання розв'язку в узагальнений степеневий ряд.

            Можна застосовувати для порівняльно простих задач (таких як лінійні)

            Чисельні методи

            це алгоритми обчислення наближений (іноді - точних) значень шуканого розв'язку y(x) на деякій вибраній сітці значень аргументу xn. Розв'язок при цьому отримується у вигляді таблиці.

            Чисельні методи не дозволяють отримати загальний розв'язок рівняння (2), вони можуть надати тільки якийсь частковий розв'язок. Вони можуть бути застосовані до широкого класу рівнянь і до всіх типів задач. Чисельні методи можна застосувати тільки до конкретно поставлених задач. В цьому сенсі задача повинна бути добре обумовлена, тобто малі зміни початкових умов не повинні приводити до великих змін в розташуванні інтегральних кривих.

            Метод розкладання в узагальнений степеневий ряд (метод Коші - Тейлора)

            Теорія - див. Демідович, Марон с128-134

            Приклад: image-20220423195533835 Розв'язок шукаємо у вигляді ряда Маклорена, (частинний випадок ряда Тейлора) image-20220423195612020 З початкових умов: image-20220423195636978

            image-20220423203204537 тобто отримали розкладення: image-20220423203224437 Якщо початкові умови задаються в довільній точці x=x0, то розв'язок шукається у вигляді ряда Тейлора image-20220423203316567 Метод рядів Тейлора є універсальним для любого диференціального рівняння з любими початковими умовами. Однак, отримати загальний вигляд члена ряда (8) неможливо (а про збігання ряда можна судити тільки по членам ряда). Ряд (8) швидко збігається в околі початкової точки. Яцщо треба отримати розв'язок на значному інтервалі, де ряд може розбігатися, або збігатися дуже повільно, цей метод не застосовується. Самостійного значення не має.

            Метод послідовних наближень Пікара

            Теорію чистити в книжці Демидович, марок стор 134-140 Це - наближений метод. Розглянемо на практиці: image-20220423203657125 (Розв'язок цього рівняння не виражається через елементарні функції) Застосуємо формулу метода Пікара image-20220423203807316 Метод Пікара доцільно застосовувати, якщо інтеграли вдається виразити через елементарні функції.

            Методи Рунте-Кутти

            Метод ломаних Ейлера

            Розглянемо на прикладі рівняння першого порядку: image-20220430212249559

            Треба скласти таблицю значень функції розв'язку y(x) в точках x1=x0+h , x2=x0+2h , x3=x0+3h,....

            Розглянемо вузол xk=x0+kh x=xk

            тоді yk=f(xk,yk)

            Зробимо заміну похідної кінцевою різницею:

            image-20220430212525763Формула (13) - це робоча формула метода Ейлера Отримаємо геометричну інтерпретацію метода. Розглянемо інтервал [xk,xk+1] на інтегральній кривій:

            image-20220430212946474

            проведемо дотичну до інтегральної кривої в точці xk: image-20220430213030439 Тобто маємо формулу метода Ейлера (13)

            Таким чином відрізок інтерполяційної кривої на інтервалі [xk,xk+1] замінюється відрізком прямої, яка співпадає з дотичною до інтегральної кривої в точці xk тобто маємо процес, зображений на рис 2. image-20220430213211038 В методі Ейлера інтегральна крива апроксимується ламаною і метод називається методом ламаних. Отримаємо методом Ейлера розв'язок задачі Коші (9) на інтервалі x[0,1] з кроком h=0,25 image-20220430213325825 Результати розрахунку представимо в табличному вигляді. В останньому стовпчику наведемо "точний" розв'язок yk, отриманий методом Пікара за формулою image-20220430213453536 image-20220430213502007 Отримали у(1,0)=0.2204. точне значення за методом Пікара y(1,0)=0,3502 значно відрізняється від отриманого. похибка складає image-20220430213615969 Покажемо, що метод Ейлера є методом 1го порядку точності. Розкладемо функцію розв'язку у(х) в ряд Тейлора в околі точки x=xk: image-20220430213714076 Порядок членів, що відкидуються, дорівнює 2. таким чином, сам метод має перший порядок.

            Формула метода ейлера для векторного рівняння (системи рівнянь)

            image-20220430213821753

            Приклад:

            image-20220430213835580

            image-20220430213844347

            Зауваження:

            Разом із значеннями функції розв'язку в вузлах xk ми маємо можливість отримувати значення похідної функції розв'язку zk+1=yk+1

            Метод Ейлера з напівкроком

            Метод розглянемо на прикладі рівняння першого порядку: image-20220430214032687 Формулу методу запишемо у вигляді image-20220430214050477 де image-20220430214059136 Таким чином, інтегральна крива апроксимується відрізками прямих, напрямок яких співпадає з напрямком дотичної до інтегральної кривої в середній точці інтервалу [xk,xk+1] Формула (15) записується двома формулами: image-20220430214226127 Застосуємо формули метода Ейлера з напівкором (16) до розв'язання image-20220430214256248 З початковою умовою у(0)=0 Отримаємо розв'язок на інтервалі [0,1] З кроком h=0,25

            Результати представлені в таблиці. На кожному кроці використовуються дві формули (16): Спочатку перша, потім друга. image-20220430214759491 Отримали: y(1,)=0,3405 Похибка результату: Δ=0,35020,3405=0,0097 значно менша похибки метода Ейлера. В подальшому ми покажемо, що метод Ейлера з напівкроком - це метод другого порядку. Формула для векторного рівняння (або системи) має вигляд: image-20220430214951443

            Модифікований метод Ейлера

            Розглянемо на прикладі рівняння першого порядку. формула метода має вигляд: image-20220430215033828 Таким чином, на інтервалі [xk,xk+1] напрямок руху співпадає з напрямком прямої з кутовим коефіцієнтом image-20220430215133824 тобто image-20220430215145314 Це - наявний метод, бо в правій частині значення image-20220430215209294 залежить від yk+1, яке ще не визначене. Формула (18) представляється у вигляді двох формул: image-20220430215309461 Перша формула (19) називається предиктор, друга - коректор. На кожному кроці застосовується обидві формули, при чому перша - один раз, друга може застосовуватися декілька разів у вигляді image-20220430215411919 Умова закінчення розрахунку на інтервалі image-20220430215432714 Де Е - задається. Буде показано, що модифікований метод Ейлера - це метод 2го порядку. Для системи рівнянь формул метода мають вигляд: image-20220430215526175

             

            Загальна схема методів Рунге-Кутти методи Рунге-Кутти 2го порядку

            Застосуємо формулу Лагранжа про кінцеві прирости на відрізку [xk,xk+1]: image-20220430215642940 Позначимо image-20220430215653982 тоді image-20220430215705706 Диференціальне рівняння image-20220430215723087 замінюється різницевим рівнянням image-20220430215745733

            Ідея Рунге-Кутти - введення в різницеву схему (3) низку додаткових параметрів, які уточнюють наближення x. Це означає, що будується апроксимація image-20220430215857135

            Представимо (3_ у вигляді

            image-20220430215943317 Тут р1,р2,...,р3,d1,d2,..,de1 - параметри, вибором яких можна забезпечити необхідний порядок апроксимації image-20220430220036728

            Наприклад, якщо у випадку бе=1 взяти р1=1, k1=f(xk,yk), будемо мати метод Ейлера (першого порядку): image-20220430220227571 Розглянемо схеми другого порядку. Покладемо в (4) е=2 Маємо: image-20220430220338176 Покажемо, що при е=2 можна так вибрати р1,р2, щоб схема (6) мала другий порядок апроксимації. З (5) для к1 і к2 маємо: image-20220430220424412 Скористаємось розкладанням Тейлора: image-20220430220501222

             

            Розкладемо функцію двох змінних f(x,y) --- image-20220430220612026 в ряд тейлора, отримуємо лінійні члени по h: image-20220430220657008 Для того, щоб схема (16) мала другий порядок по h, необхідно, щоб одночасно коефіцієнти при yk і при yk дорівнювали 0 image-20220430220811302

             

            Тобто image-20220510151630232

            Всі схеми при image-20220510151650977, для яких має місце умова (17), мають другий порядок точності.

            Розглянемо випадок, коли p2=1;α1=1/2;p1=0 Маємо: image-20220510151811326

            Це розглянутий раніше метод Ейлера з напівкроком.

            Тепер візьмемо image-20220510151904937 Маємо: image-20220510151918125Це розглянутий раніше модифікований метод Ейлера.

            Таким чином ми показали, що метод Ейлера з напівкроком і модифікований метод Ейлера - це методи Рунге-Кутти 2го порядку.

            Методи Рунге-Кутти 3го і 4го порядків

            Методи 3го порядку

            Загальна схема методів 3го порядку: image-20220510152057150 Оскільки на момент обчислення к3 вже відомі к1 і к2, схему можна деталізувати: image-20220510152131467 Як і у випадку методів 2го порядку, можна показати, що калетним(?) вибором image-20220510152222491 можна забезпечити схемі (1)-(2) третій порядок точності.

            Умова цього наступна: image-20220510152301762

            Серед безліч схем 3го порядку наведемо наступні

             

            Запишемо формули методів А)-Г) для векторного рівняння:

            А)image-20220510152530326

            Б) image-20220510152547599

            В) image-20220510152556809

            Г) image-20220510152608507

            Приклад

            Розглянемо систему: image-20220510152633014 Запишемо формули методу А): image-20220510152653325 image-20220510152708351

            Методи Рунге-Кутти 4го порядку

            Цей метод найбільш часто використовується в технічних задачах. Схема метода для рівняння першого порядку: image-20220510152820589 Значення image-20220510152828492 отримані з умови співпадання різницевої схеми (8) з рядом тейлора включно до членів порядку h4

            Для цієї схеми image-20220510152935821 Схема метода Рунге-Кутти 4го порядку для векторного рівняння image-20220510153003669:

            image-20220510153017164

            Приклад

            Система: image-20220510153034769

            image-20220510153048546

            Висновки по методам Рунге-Кутти

            1. Всі методи є явними (використовують значення аргументу x[xk,xk+1))
            2. Методи є однокроковими (використовують інформацію про у і f(x,y) тільки з інтервалу [xk,xk+1)
            3. Припускають обчислення із змінним кроком
            4. легко переносяться на системи диференціальних рівнянь,

            Підвищення точності чисельних методів інтегрування звичайних диференціальних рівнянь Правило Рунге

            Розглянемо диф рівняння image-20220510161125869 Нехай відомі значення точного розв'язку image-20220510161146993 а наступне значення image-20220510161210331 розкладене по степеням кроку h в ряд Тейлора: image-20220510161242093

            Застосуємо чисельний (наближенний) метод до точних значень image-20220510161314118 і отримаємо image-20220510161325297 Розкладемо його в ряд Тейлора: image-20220510161341590

            Якщо для всіх n виконується рівність image-20220510161412853 причому image-20220510161420823 для деякого n, то число m називається порядком точності наближеного методу.

            При цьому похибка методу на кожному кроці має порядок O(hm+1). Якщо число m (порядок) відомо, то для грубої оцінки похибки методу можна застосувати принцип Рунге

            Розрахунок проводять на двох сітках значень аргументу х. Для першої сітки відстань між вузлами становить h, для другої - 2h. Позначимо точний розв'язок рівняння (1), як y(x). Наближений розв'язок отриманий на густій сітці з кроком h позначимо, як image-20220510161833413 . Тоді на одному кроці обчислень маємо оцінку:

            image-20220510161900284

            вважаючи А постійним для кожного кроку обчислень, зробимо 2n кроків по густій сітці. Отримаємо: image-20220510161944667

            Тепер зробимо n кроків по розрідженій сітці (H=2h). Отриманий наближений розв'язок позначимо, як image-20220510162026558

            image-20220510162033218

            Для цього розв'язку маємо: image-20220510162103108

            Формула (6) дозволяє уточнити розв'язок, отриманий на густій сітці за рахунок розв'язку, що отриманий на розрідженій сітці. Враховуючи (4), запишемо уточнений розв'язок image-20220510162649972 Підставимо 6 в 7: image-20220510162704528 Формула (8) уточнює розв'язок image-20220510162724868 в вузлах, що співпадають, поправка складає image-20220510162751853

            Значення поправок в решті вузлів густої сітки знайдемо на основі найпростішої інтерполяції: image-20220510162831737

            Приклад

            image-20220510162846588

            10.05.22

            Багатокрокові методи розв'язання задачі Коші. Метод Адамса 4го порядку

            В раніше розглянутих методах Рунге-Кутти значення yk+1 обчислювалося на основі апрокимацій функції f(x,y) де x[xk,xk=1]. тобто використовувалися х і у, що відносилися тільки до одного кроку. Такі методи називаються однокроковими. Логічно припустити, що більшу точність можна отримати, якщо використати інформацію з попередніх кроків:

            Саме на цій ідеї основані багатокрокові мтеоди.

            Розглянемо лінійні багатокрокові методи.

            Загальна формула

            Формула 1 дає опис загального лінійного р-крокового методу. Якщо β0=0 - метод називається явним, інакше маємо неявну схему. В цьому випадку для обчислення yk+1 треба знати fk+1, яке невідоме.

            При застосуванні багатокрокових методів на кожному кроці обчислень звичайно використовують одночасно дві схеми - явну і неявну. Явна схема називається предиктор, неявна - коректор. Явна схема дає початкове наближення yk+10, неявна - уточнює це наближення.

            Метод Адамса (Штермера)

            Метод Адамса отриманий за допомогою другої інтерполяційної формули Ньютона, в якій утримуються кінцеві різниці до 4го порядку включно

            Робочі формули методу

            1. Явна схема:
            2. Неявна схема:

            Схема рішення:

            1. використати предиктор 2 для обчислення початкового наближення yk+10
            2. покласти i=0. По формулі коректора 3 знайти уточнення

            Зауваження

            1. Щоб розпочати обчислення за формулою 2 необхідно мати значення Для їх знаходження треба використати інший метод 4го порядку, наприклад, метод Рунге-Кутти. Тобто спочатку треба тричі використати метод Рунге-Кутти 4го порядку з кроком h, а тільки потім формулу 2. Така процедура називається розгонною

            2. Схему коректора можна використовувати на кожному кроці декілька разів. За таких умов кінець розрахунків на кожному кроці задається нерівністю:

               

            У випадку, коли маємо систему диференціальних рівнянь, формули Адамса набувають вигляду:

            Запишемо розрахункові формули метода Адамса для системи:

            Предиктор:

            Коректор:

            Формула 6 (? про систему чи шо) дозволяє уточнити розв'язок, отриманий на густій сітці за рахунок розв'язку, що отриманий на розріженій сітці. Враховуючи 4, запишемо уточнений розв'язок

            Підставимо 6 в 7:

            Формула 8 уточнює розв'язок y~(x,h) в вузлах, що співпадають, поправка складає

            Значення поправок в решті вузлів густої сітки знайдемо на основі найпростішої інтерполяції:

            Приклад:

            В таблиці останній стовпчик y(x) - точний розв'язок, отриманий раніше методом Пікара.

            Використаний метод чисельного інтегрування - метод Ейлера (m=1), тобто формула 9 набуває вигляду:

            Другий стовпчик - наближений розв'язок з кроком h=0.25

            Третій стовпчик - наближений розв'язок з кроком 2h=0.5

            Четвертий стовпчик: добавка до розв'язку: в співпадаючих вузлах вираховується за формулою 10, в неспівпадаючих вузлах - за формулою 9.1

            П'ятий стовпчик - уточнений розв'язок

            Висновок

            Уточнений розв'язок точніший за розв'язок y~(x,h)

            Жорсткі рівняння. Неявний метод Ейлера

            Попередні відомості

            Постійна часу в розв'язку диференціального рівняння - це час, необхідний для його зменшення в 12 разів.

            Приклад

            Маємо рівняння

            Отримаємо загальний розв'язок:

            Стійке диференціальне рівняння називається жорстким, якщо воно має частинний розв'язок у вигляді експоненти, що складає, постійна часу якого дуже мала у зрівнянні з довжиною інтервала, на якому шукається цей розв'язок. Для системи y¯=f¯(t,y¯)

            Швидкості спадання пов'язані з власними значеннями матриці якості

            Розглянемо наспупний приклад: Знайдемо розв'язок системи

            Матриця Якості маж вигляд:

            Знайдемо власні значення:

            Нескладно бачити, що λ1=1;λ2=1000 Запишемо розв'язок системи 3 у вигляді:

            Знайдемо . Для цього продиверенціюємо 4:

            З початкових умов маємо:

            з 3 і 5 при t=0 отримаємо:

            Таким чином

            Отже розв'язок має вигляд:

            Знайдемо ектстремальні значення розв'язку 6

            Спробуємо знайти чисельний розв'язок системи 3 методом Ейлера: Маємо робочі формули:

            Нехай h=0.01 зробимо один крок

            Отримали, що вже після першого кроку

            Розглянемо неявний метод Ейлера:

            Який для системи 3 записуються так:

            Формули 7 приводяться до вигляду:

            Зробимо один крок за схемою 8

            Система 9 дає

            Ці значення більш прийнятні, ніж отримані звичайним методом Ейлера Основне протиріччя жорстких систем полягає в тому, що малий крок інтегрування, потрібний для відтворення щвидкоплинних процесів не може бути збільшений в подальшому, коли похідна розв'язку стає суттєво меншою за її початкове значення. Перший участок (на рисунку) з швидким зростанням (спаданням) розв'язку називається пограничним шаром. Для розглянутої системи це t[0;0.01]

            Властивість жорсткості

            Задача Коші

            Називається жорсткою в інтервалі t[a,b], якщо відношення

            Де λi - власні значення матриці Якості , в яку підставлено розв'язок y(t) в точці t

            Величину μ(t) можна назвати коефіцієнтом жорсткості системи 10 Для розглянутої систему 3 μ=10001>>1, тобто система 3 - жорстка

            17.05.22

            Чисельне методи розв'язання крайових задач

            Загальні відомості про крайові задачі

            Постановка задачі:

            Для заданого диференціального рівняння

            треба знайти розв'язок y=y(x),x[a,b], який задовольняє крайовим умовам:

            Крайові задачі поділяють на лінійні і нелінійні. Лінійні задачі мають місце, коли диференціальне рівняння лінійне, і лінійні крайові умови. Лінійні задачі поділяються на однорідні і неоднорідні

            Приклад

            Найпростіша двоточкова крайова задача. Знайти функцію y=y(x), яка задовольняє диференціальному рівнянню другого порядку

            І приймає при x=a і x=b(a<b) задані значення y(a)=A;y(b)=B

            Геометрично це означає, що треба знайти інтегральну криву рівняння 3, яка проходить через дані точки M(a,A) і N(b,B)

            Крайова задача може:

            а) не мати розв'язків Б) мати єдиний В) мати декілька, або нескінчено багато

            Приклад 2

            Диференціальне рівняння

            (4)y+y=0

            має загальний розв'язок

            А) нехай y(0)=0;y(π)=0 тоді c1=0,c2 - довільне, тобто маємо нескінчену множину розв'язків y=csin(x)

            Б) Нехай y(0)=1;y(π)=1 тоді отримаємо тобто розв'язків задача не має

            В) Нехай y(0)=1;y(π2)=1, тоді і маємо єдиний розв'язок

            Редукція (зведення) лінійної крайової задачі до задачі з початковими умовами

            Розглянемо метод на прикладі рівняння другого порядку

            З крайовими умовами:

            Спочатку розглянемо однорідне рівняння:

            1. задамо початкові умови Розв'язуючи чисельно задачі Коші 7 8 знайдемо розв'язок u1(x)

            2. Задамо інші початкові умови

               

            і знайдемо розв'язок задачі 7 9 u2(x) Можна показати, що розв'язки u1(x) і u2(x) лінійно незалежні

            Тоді маємо загальний розв'язок рівняння 7 у вигляді

            Тепер знайдемо розв'язок рівняння 5

            З довільними початковими умовами, наприклад

            і знайдений розв'язок рівняння 5 буде мати вигляд:

            Залишається визначити с1 і с2 з крайових умов 6 З 13 початкових умов 8 і 9 маємо де u1(b),u2(b),v(b) знайдені в результаті рішення відповідних задач Коші Далі диференціюємо 13

            Звідки маємо:

            Значення u1(b),u2(b),v(b)також знайдені в результаті розв'язування задач Коші. Таким чином отримали як залежності від с1 с2. Підстановка цих залежностей в крайові умови 6 призводить до системи лінійних рівнянь відносно с1 і с2:

            Розв'язуючи яку знайдемо невідомі с1 і с2 таким чином крайова задача звелась до розв'язання трьох задач Коші і системи лінійних рівнянь 17

            Зведення задачі до рівнянь в кінцевих різницях

            Розглянемо рівняння другого порядку

            Де права частина в загальному випадку - нелінійна функція, і крайові умови 6 Розіб'ємо інтервал [a,b] точками x0=a, x1,x2,..,xn1,xn=b на n рівних відрізків довжиною . Треба знайти

            Запишемо ряди Тейлора

            Віднімемо 20 від 19, маємо

            Звідкіля отримаємо:

            Додаючи 19 до 20 маємо:

            Таким чином для x=xk з 18 21 22 маємо з точністю до 0(h2)

            В системі 23 всього n-1 рівняння відносно невідомих y0,y1,..yn, яких n+1 Треба ці дві умови знайти. Їх знайдемо з граничних умов 6 Для цього запишемо ряди Тейлора:

            Домножимо 24 на 4 і віднімемо від нього 25:

            Для правого кінця також маємо:

            Отримаємо з 27 28:

            Підставимо вирази 26 29 в крайові умови і отримаємо з врахуванням 23 наступну систему рівнянь:

            Таким чином задача 18, 6 звелася до системи n+1 рівнянь з n+1 невідомими. розглянемо випадок лінійного рівняння

            Кінцево-різницевий аналог цього рівняння

            Тоді система 30 приймає вигляд:

            Це - трьохдіагональна система лінійних рівнянь. Її можна розв'язати звичайним методом Гайса, Але більш економічним методом її розв'язування є розглянутий раніше метод прогонки, який застосовується при інтерполяції сплайнами.

             

            Метод прогонки

            При застосуванні методу кінцевих різниць до крайових задач для диференціальних рівнянь другого порядку маємо трьохдіагональну систему лінійних алгебраїчних рівнянь, кожне з яких містить три сусідні невідомі. Для розв'язування такої системи розроблений спеціальний метод - метод прогонки.

            Розглянемо лінійне диференціальне рівняння

            З двоточечними лінійними крайовими умовами

            - неперервні функції на [a,b]

            Як було показано в попередній лекції, замість диференціального рівняння 1 отримаємо систему кінцево різницевих рівнянь:

            Після нескладних перетворень запишемо 3 у вигляді

            Де позначено:

            Для похідних на кінцях x0=a;xn=bберемо односторнні похідні

            Звідсіля згідно умов 2 отримаємо

            Розглянемо метод прогонки: Розв'яжемо систему 4 відносно yi маємо:

            Нехай за допомогою повної системи 4 6 з рівняння 7 виключемо невідоме yi1' тоді це рівняння має вигляд

            де - деякі коефіцієнти. 3 8 маємо:

            Підставляючи цей вираз в рівняння 4 отримаємо:

            Порівнюючи 8 і 9 отримаємо для визначення с1 і d1 рекурентні формули

            Визначимо окремо c0,d0 з першої крайової умови 6 отримаємо

            З іншого боку з формули 8 при i=0 маємо

            Порівнюючи два останні рівняння знаходимо:

            На основі формул 10 11 послідовно визначаються коефіцієнти

            (прямим ходом)

            Зворотний хід

            починається з визначення yn Використовуючи другу крайову умову 6 і формулу 8 при i=n-1 отримаємо систему двох рівнянь

            Розв'язуючу її відносно yn, будемо мати

            Тепер за формулою 8 послідовно знайдемо

            Для найпростіших крайових умов

            Формулм для c0,d0,y0,yn спрощуються Дійсно, якщо з формул 12 маємо:

            Треба відмітити, що метод прогонки має стійкий обчислюваний алгоритм і похибки округлень не викликають нескінченного зростання похибки рішення

            Приклад

            Методом прогонки знайти розв'язок крайової задачі

            Рішення

            Візьмемо h=0.1, перейдемо від рівняння 17 до кінцеворізницевого аналого:

            Згідно формули 15 маємо

            24.05.22

            Лінійна прогонка із застосуванням спеціальної функції

            Розглянемо лінійне диференціальне рівняння

            (1)y+p(x)y+q(x)y=f(x)

             

            З найпростішими крайовими умовами

            (2)y(a)=A;y(b)=B

             

            Введемо спеціальну функцію вигляду:

            (3)y=τ(x)y+S(x)

            Де τ(x);s(x) - деякі неперервні функції на [a,b] Нехай τ(x);s(x) також неперервно диференційовані на [a,b] Продиференціюємо 3, отримаємо

            (4)y=ry+ry+sr=r(x),r=r(x),s=s(x)

            Підставимо замість y вираз 3:

            (5)y=ry+r(ry+s)+s=y(r+r2)+rs+s(s)

             

            Рівність 5 означає, що y(x)​ розкладено в рід по степеням розв'язку y(x) З 1 маємо:

            (6)y=pyqy+f

            Замінемо в 6 y виразом 3:

            (7)y=p(ry+s)qy+f==y(prq)ps+f

             

            Порівняємо розкладення 5 і 7 Прирівнюючи коефіцієнти при однакових степенях у маємо систему:

            (8){rs+s=ps+fr+r2=prq

             

            Тобто отримали систему двох диференціальних рівнянь першого порядку відносно τ(x)​ і s(x)​:

            (9){s+s(r+p)f=0r+r(r+p)+q=0

             

            З крайових умов 2 отримаємо початкові умови для інтегрування системи 9 для x=b з 3 отримаємо:

            y(b)=r(b)y(b)+S(b).

            Оскільки y(b) невідоме, покладемо

            (10)r(b)=0

            Тоді:

            (11)S(b)=y(B)=B

             

            Таким чином отримали початкові умови для інтегрування системи 9 в зворотному часі

            (12)r(b)=0;S(b)=B

            Інтегруючи систему 9 від x=b​ до x=a​ з кроком h вибраним чисельним методом, отримаємо значення τ(a),s(a) Тепер запишемо початкове рівняння 1 в формі Коші

            (13){y=zz=pzqy+f

             

            З виразу 3 отримаємо другу початкову умову для інтегрування системи 13

            y(a)=r(a)y(a)+s(a)==r(a)A+s(a).

            Виконаємо інтегрування системи 13 з початковими умовами

            y(a)=A,y(a)=z(a)=Ar(a)+S(a)

            Від x=a до x=b і тим самим розв'яжемо початкову крайову задачу. Таким чином крайова задача звелася до двох задач з початковими умовами (задач Коші). При цьому інтегрування системи 9 в зворотному напрямку можна інтерпретувати, як зворотню прогонку, а інтегрування системи 13 - як пряму прогонку

            P.S

            В деяких випадках процес інтегрування системи 9 може бути нестійким. тоді можна пропонувати інший зв'язок між y(x) і y(x), а саме:

            Розглянемо процес лінійної прогонки для цього випадку. Диференціюєсо 14:

            Розв'яжемо 15 відносно старшої похідної

            Скористаємося виразом 6, підставивши в нього замість у вираз 14

            Вирази 16 і 17 - це розкладання функції y(x) в ряд по степеням y(x) прирівнюючи коефіцієнти при однакових степенях, маємо:

            Отже, маємо систему диференціальних рівнянь

            Отримаємо початкові умови для інтегрування системи 18 з урахуванням вигляду крайових умов 2 і спеціальної функції 14:

            Оскільки y(a) - невідоме, покладемо коефіцієнт при y(a)

            Отже інтегруємо систему 13 з початковими умовами 20 21 від x=1 до x=b з кроком h, і в результаті знайдемо τ(b) і s(b) Тепер можна знайти другу початкову умови для інтегрування системи 13:

            Інтегрування системи 13 здійснюємо в зворотному часі від x=b до x=a з кроком h і отримаємо розв'язок початкової крайової задачі.

            Приклад застосування методу прогонки з використанням спеціальної функції

            Знайти розв'язок диференціального рівняння З крайовими умовами На інтервалі x[0,12] з кроком 0.1 Введемо спеціальну функцію

            Оскільки y(1/2) невідоме, то покладемо τ(1/2)=0, тоді з 2 отримаємо:

            З 3 маємо:

            Записуємо диференціальні рівняння відносно τ(x) і s(x). Вони мають вигляд:

            Початкові умови для інтегрування системи 5:

            Використаємо для інтегрування метод Ейлера. Результати представлені в табл 1. Інтегрування відбувається в зворотньому часі від х=0.5 до х=0 з кроком h=-0.1

            Табл1

            Робочі формули метода Ейлера для системи 5:

            В таблиці позначено:

            Після інтегрування отримали:

            Тепер з 2 можна визначити

            Це - друга початкова умова для інтегрування рівняння 1 від х=0 до х=0.5 з кроком h=0.1. Перша початкова умова y(0)=1 Інтегруємо рівняння 1 у вигляді системи:

            Використаємо метод Ейлера. Результати представлені в табл. Робочі формули метода Ейлера для системи 7

            Позначимо:

            Табл 2

            Zn=z(0)=y(0)=4.257

            Отже отримали розв'язок крайової задачі який представлений другим та третім стовпчиками табл2

            Перевірка розв'язку

            Крайова умова y(0.5)=0.5 не виконана, бо отримали y(0.5)=0.019. Похибка обумовлена вибраним методом першого порядку і великим кроком чисельного інтегрування.

            Метод стрільби (балістичний) для систем диференціальних рівнянь

            Ще один спосіб зведення крайової задачі до задач Коші дає так званий балістичний метод, або метод стрільби.

            Розглянемо крайову задачу для системи двох диференціальних рівнянь першого порядку з нелінійними крайовими умовами:

            (5.1)u(x)=f1(x,u,v),v(x)=f2(x,u,v)x(a,b)
            (5.2)φ[u(a),v(a)]=0,ψ[u(b),v(b)]=0

            Виберемо довільне значення u(a)=η і розглядатимемо першу з крайових умов (5.2) як алгебраїчне рівняння φ[η,v(a)]=0 звідки визначимо v(a)=ζ(η) Візьмемо значення u(a)=η,v(a)=ξ, за початкові умови задачі Коші для системи (5.1). Зінтегрувавши цю задачу, одержимо розв’язок, залежний від параметра: u=u(x,η),v=v(x,η) Знайдений розв’язок справджує першу з крайових умов (5.2), але, взагалі кажучи, не справджуватиме другу з цих умов.

            Отже, параметр η слід підібрати таким чином, щоб

            (5.3)ψ¯(η)ψ[u(b,η),v(b,η)]=0

            Знайти корінь алгебраїчного рівняння (5.3) можна різними наближеними методами.

            Метод дихотомії.

            Робимо «пробні постріли» – розрахунки з довільним чином вибраними значеннями ηi допоки серед величин ψ(ηi) не буде різних за знаком. Пара відповідних їм значень параметра ηi та ηi+1 утворять першу «вилку», в межах якої міститься корінь рівняння (5.3). Поділяючи «вилку» послідовно навпіл до одержання потрібної точності, здійснюємо «пристрілку» параметра η. Саме завдяки цьому процесу метод і одержав назву методу стрільби.

            Однак знаходження кожного нового значення функції ψ¯(η) вимагає числового інтегрування системи (5.1), тобто є досить громіздким. Тому для відшукання кореня рівняння (5.3) доцільно застосовувати більш швидкі методи.

            Метод січних.

            Перші два «постріли» η0 та η1 робимо навмання, а наступні значення обчислюємо згідно з формулою січних:

            (5.4)ηi+1=ηi(ηiηi1)ψ¯(ηi)ψ¯(ηi)ψ¯(ηi1),i=1,2,

            Можна використовувати й інші подібні методи.

            Особливо просто за допомогою методу стрільби розв’язуються лінійні крайові задачі. Справді, нехай задача (5.1),(5.2) є лінійною, тобто має вигляд

            (5.5)u(x)=p1(x)u+q1(x)v+r1(x)v(x)=p2(x)u+q2(x)v+r2(x),x(a,b)

            (5.6)α1u(a)+β1v(a)=γ1,α2u(a)+β2v(a)=γ2

            Тоді початковими значеннями відповідної задачі Коші будуть

            (5.7)u(a)=η,v(a)=ξβ11(γ1α1η)

            Легко бачити, що розв’язок задачі Коші (5.5),(5.7) лінійно залежатиме від параметра η, тому функція ψ¯(η) також буде лінійною. Але лінійна функція одного аргументу цілком визначається своїми значеннями в довільних двох точках η0 і η1 , а її графіком є пряма, тобто вона співпадає зі своєю січною. Тому знайдене за формулою (5.4) значення η2 є точним коренем рівняння (5.3), тобто розрахунок із цим значенням дає шуканий розв’язок крайової задачі (5.5),(5.6). Отже, для розв’язання лінійної крайової задачі достатньо трикратного інтегрування задач Коші.

            Зауважимо, що для лінійних задач можна суттєво зменшити обсяг роз-рахунків, врахувавши структуру загального розв’язку системи (5.5). Як відомо з курсу звичайних диференціальних рівнянь, загальний розв’язок лінійної неод-норідної системи рівний сумі загального розв’язку відповідної однорідної системи та деякого частинного розв’язку неоднорідної системи.

            Знайдемо частинний розв’язок системи (5.5), що відповідає значенню η0=0 і позначимо його через u0(x),v0(x) Розглянемо задачу Коші для відповідної однорідної системи вигляду

            u(x)=p1(x)u+q1(x)v,x(a,b);v(x)=p2(x)u+q2(x)v,u(a)=η1=1,v(a)=α1β11

            знайдемо її розв’язок і позначимо його через u1(x),v1(x) Тоді загальним розв’язком неоднорідної системи, який справджуватиме першу з крайових умов (5.6), буде однопараметрична сім’я вигляду

            (5.8)u(x)=u0(x)+Cu1(x),v(x)=v0(x)+Cv1(x)

            Значення параметра С знаходимо з другої з крайових умов (5.6):

            C=α2u0(b)+β2v0(b)γ2α2u1(b)+β2v1(b).

            Підклавши знайдене значення С в (5.8), одержимо розв’язок крайової задачі (5.5),(5.6).

            Метод стрільби є простим і успішно застосовним до більшості крайових задач типу (5.1),(5.2), а також до звичайних диференціальних рівнянь: адже, як відомо, будь-яке диференціальне рівняння п-го порядку можна звести до системи п рівнянь першого порядку. Труднощі виникають тоді, коли крайова задача (5.1),(5.2) добре обумовлена, а відповідна їй задача Коші погано обумов-лена. При цьому чисельне інтегрування задачі Коші визначає функцію ψ¯(η) з великою похибкою, що ускладнює організацію ітерацій. У такому випадку пробують задавати початкові умови на правому кінці заданого проміжку x=b і інтегрувати справа наліво; нерідко при цьому стійкість схем покращується. Якщо й це не допомагає, то таку крайову задачу слід розв’язувати або спеціаль-ними, або різницевими методами.

            ПРИКЛАД

            За допомогою методу стрільби знайти розв’язок крайової задачі:

            (5.9)y2x2y=2x1,x(1;2)
            (5.10)9y(1)5y(1)=0,y(2)y(2)+1=0

            Розв’язання

            Задача (5.9),(5.10) є лінійною, тому згідно з алгоритмом методу стрільби для відшукання її розв’язку достатньо інтегрування трьох задач Коші для рівняння (5.9) з початковими умовами y(1)=ηi,y91)=1.8ηi,i=0,2,.. ; тут значення похідної в точці x=1 визначені з першої з крайових умов (5.10). Перші два «постріли» η0 та η1 робимо навмання, беручи довільні два числа, а значення η2 обчислюємо за формулою січних (5.4), в якій з урахуванням другої з крайових умов (5.10) слід узяти ψ¯(ηi)=yi(2)yi(2)=1 де yi(x) – розв’язок задачі Коші для рівняння (5.9) з параметром ηi

            Нехай зроблені «постріли», наприклад, η0=5,η1=0 Зінтегрувавши рівняння (5.9) спершу при початкових умовах y0(1)=5,y0(1)=9 і вдруге – при початкових умовахy1(1)=0,y1(1)=0 одержимо розв’язки відповідних задач Коші:

            y0(x)=23(8x2+x1)x,y1(x)=13(2x2+x1)x

            Тоді згідно з формулою січних (5.4)

            (5.11)η2=η1(η1η0)ψ¯(η1)ψ¯(η1)ψ¯(η0)=5(y1(2)y1(2)+1)y1(2)y1(2)+1(y0(2)y0(2)+1)

            Підклавши в (5.11) знайдені функції y0(x) та y1(x) одержимо значення η2=5 Зінтегрувавши рівняння (5.9) втретє, цього разу при початкових умовах y2(1)=5,y2(1)=9 одержимо розв’язок відповідної задачі Коші

            (5.12)y2(x)=x(4x+1)

            який згідно з алгоритмом методу стрільби для лінійних крайових задач повинен бути розв’язком крайової задачі (5.9),(5.10). Безпосереднім підкладанням легко переконатися, що функція (5.12) дійсно справджує обидві крайові умови з (5.10). Отже, розв’язок крайової задачі (5.9),(5.10) y(x)y2(x)=x(4x+1)